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INTRODUCTION 


It  is  well  known  that  interatomic  interactions  play  a  crucial  role  in 
the  determination  of  macroscopic  phenomena  such  as  plastic  flow,  fracture  and 
diffusion-controlled  processes  like  void  or  precipitate  formation.  In  many 
cases,  metallurgical  properties  of  practical  concern  are  controlled  by  inter¬ 
actions  involving  defects,  often  point  defects  such  as  vacancies  or  impurity 
ions.  A  theory  of  the  properties  of  defects  in  terms  of  interatomic  inter¬ 
actions  is  therefore  vital  to  a  fundamental  understanding  of  important  metal¬ 
lurgical  processes. 

Almost  all  studies  of  the  structures  and  properties  of  defects  in  metals 
are  based  on  the  assumption  that  interatomic  forces  are  derived  from  central, 
pairwise  interatomic  potentials.  This  assumption  specifically  excludes  the 
possibility  of  multi-ion  interactions,  i.e.,  interactions  that  involve  three 
or  more  ions  and  that  are  not  simply  a  sum  of  pairwise  interactions.  Although 
theory  shows  that  multi-ion  interactions  do  indeed  exist,  calculations  of  the 
interaction  energies  are  difficult,  and,  as  a  result,  not  much  was  known  about 
their  importance  until  recently.  Now,  however,  there  is  a  growing  body  of 
evidence  that  tells  us  that  three-ion  interactions  can  be  quite  strong,  and 
that  ignoring  them  can  lead  to  spurious,  sometimes  unphysical  results.  This 
is  particularly  true  of  certain  calculations  involving  defects,  where  recent 
results  indicate  that  one  must  include  at  least  three-ion  interactions,  in 
addition  to  pairwise  interactions,  to  obtain  accurate  predictions  of  some 
properties. 

Accordingly,  the  objectives  of  the  work  reported  here  were  to  (1)  assess 
the  relative  importance  of  two-  and  three-ion  interactions  involving  defects  in 
metals,  (2)  develop  and  test  methods  for  calculating  three-ion  interaction 
energies,  and  (3)  develop  simplified,  approximate  models  of  multi-ion 


interactions  for  use  in  computer  simulation  studies  of  defects.  Emphasis  was 


placed  on  calculations  of  the  properties  of  point  defects  in  aluminum  because 
three-ion  interactions  are  known  to  play  an  important  role  in  the  determina¬ 
tion  of  some  point  defect  properties,  and  because  ample  experimental  data  are 
available  for  testing  the  results  of  the  calculations. 

Section  II  of  this  report  contains  a  brief  review  of  the  basic  ideas  in 
multi-ion  interaction  theory  and  a  summary  of  progress  toward  meeting  our  ob¬ 
jectives.  Additional  technical  information  is  reported  in  reprints  and  pre¬ 
prints,  which  are  included  as  Appendices. 

II.  TECHNICAL  DISCUSSION 
A.  Background 

This  section  contains  a  brief  review  of  pseudopotential/perturbation 
theory  as  applied  to  simple  nontransition  metals.  1116  principal  purpose  here 
is  to  establish  terminology  and  notation;  comprehensive  reviews  are  available 
elsewhere  (Ref.  1,  2).* 

We  visualize  a  system  of  N  ions  with  fixed  core  states  embedded  in 
a  valence  election  gas  of  density  pQ.  The  average  energy  per  ion  can  be  ex¬ 
pressed  as  the  sum  of  a  direct  ion-ion  interaction  energy  E^  and  an  electronic 
energy  Eq.  Thus 

E  =  Ei  +  Eg  (1) 

The  electronic  energy  is  expressed  as  a  perturbation  expansion  in 
the  electron-ion  pseudopotential  W 

Ee’*Eo+Ei  +  E2  +  E3+...  (2) 

where  Eh  contains  terms  of  n’th  order  in  W.  The  crystal  pseudopotential  is 
the  sum  of  two  terms 

W  -  "ion  +  **v 

* References  cited  in  the  text  are  listed  at  the  end  of  this  report. 


(3) 


where  w^on  is,  in  general,  a  nonlocal  potential  (an  operator),  while  Wy  is  the 
self-consistent  local  potential  energy  of  a  valence  electron  in  the  presence 
of  all  other  valence  electrons.  Self-consistency  means  that  the  plane  wave 
matrix  elements  of  Wv  satisfy  the  modified  Poisson  equation  (in  ftydbergs) 

Wv(q)  =  ~7[l“G(q) ]pv(q)  (4) 

■* 

where  pv(q)  is  the  Fourier  transform  of  the  valence  electron  density  and  G(q) 
is  a  correction  for  exchange  and  correlation  interactions  among  electrons. 

The  electron  density  may  also  be  written  as  a  perturbation  eiqpansion 

Pv(q)  =  &qoPo  +  +  +  •••  (5) 

where  the  linear  term  pi(q)  is  a  sum  of  two  terms,  one  containing  a  matrix 
element  of  Wion  and  the  other  proportional  to  WyCq) .  Hie  second  order  term 
involves  an  integral  over  products  of  these  matrix  elements,  the  third  order 
term  involves  double  integrals,  etc.  To  calculate  the  self  consistent  pseudo¬ 
potential  to  order  n,  one  must  truncate  (5)  at  the  n'th  term,  substitute  in 
(4)  and  solve  for  Wytq);  for  n  >  2  this  means  solving  an  integral  equation. 

It  is  easily  shown  that,  to  calculate  Eq  to  order  n,  one  needs  Wy(q) 
to  order  n-1  (Ref.  1,  2).  A  second  order  calculation  of  the  energy  therefore 
involves  a  linear  screening  approximation  to  the  self  consistent  solution  of 
(4)  and  (5).  Hie  well  known  result  is 

<£+q|W|k>  ~  <k+q|Wion|k>/e(q)  (6) 

where  e,  the  dielectric  function,  depends  on  the  average  electron  density  as 
well  as  on  the  momentum  change  q. 


After  rearranging  terms,  the  average  energy  per  ion,  to  second  order 


in  W,  can  be  written 

B2  "  V°  +  2^  ll  Vij(Rij>  <7> 

i  i*j 

where  V0  is  independent  of  the  arrangement  of  ions,  Rf j  is  the  distance  between 
ions  at  R^  and  Rj,  and  v^j  is  the  effective  interaction  between  ions,  or  the 
interatomic  potential.  We  include  the  subscript  ij  to  account  for  the  possi¬ 
bility  that  different  types  of  ions  may  be  present.  For  a  homogeneous  system, 
the  function  V(R)  is  the  same  for  all  ion  pairs  and  can  be  expressed  as  the 
direct  ion-ion  potential  plus  the  Fourier  transform  of  a  function,  Fn(q), 
called  the  energy  wavenumber  characteristic  of  the  material. 

In  a  higher  order  approximation,  one  does  not  obtain  a  central,  pair¬ 
wise  interatomic  potential.  The  energy  can  be  written  in  a  form  similar  to 
(7),  (Ref.  3,4) 

®  =  VQ  +  1-  H  VijtR^Rj)  (8) 

*  i  i*j 

but  now  the  effective  potential  is  a  pairwise  term  plus  a  sum  of  multi-ion 
interactions.  If  all  ions  are  identical  (Ref.  3,4)  we  have 


Vij(Ri*Rj)  -  V(Ri;))  +  j  l  VstR^Rj.Rfc)  +  ••• 


The  second  term  is  a  three-ion  interaction  energy  that  results  from  carrying 
the  calculation  of  B  to  third  order;  the  next  term  would  be  a  fourth  order. 


four  ion  term,  etc 


It  is  possible  to  develop  formal  expressions  for  higher  order  ener¬ 
gies  in  terms  of  matrix  elements  of  nonlocal  pseudopotentials,  and  thus  deter¬ 
mine  the  corresponding  formulas  for  V3,  V4,  etc.  However,  even  in  third  order, 
the  result  is  too  complicated  to  be  of  use  in  practical  computations  (Ref.  5). 
It  is  therefore  necessary  to  make  local  approximation  in  third  and  higher  order 
terms,  which  amounts  to  taking 

<k+q|W|k>  =  NS(q)<k+q|u>|k>  (10) 

~  S(q)u)(q) 

-► 

Here  S(q)  is  the  structure  factor, 

s<?>  -  5  1  "» 

i 

u)  is  the  nonlocal  pseudopotential  operator  for  a  single  ion  (all  ions  are 
assumed  identical),  and  w(q)  is  a  local  approximation  to  N<k+q| w|k>.  This 
leads  to  the  following  expression  for  the  third  order  energy  (Ref.  6,7). 

e3  *  1  g(q-|,<l2'cI3)<,J(cIl)a,({I2)<jJ(<l3)s(<Jl)s(<l2)S(q3)  (12) 

*1  «3 

where  5^,  is  the  average  volume  per  ion,  q3  =  ”(qi+q3)»  and  g(qvq2'33)  is  the 
function  defined  by  Lloyd  and  Sholl  (Ref.  6).  If,  for  convenience,  we  also  use 
the  local  approximation  in  the  second  order  term  (it  is  not  necessary  to  make 
this  approximation)  then  the  effective  potential  in  (8)  becomes 


q 

■>  + 

+  1?°  I  w(qi)e“iqi  Ri3  £  g(qi,q2/q3)w(q2)w(q3) 


+  “2“  l  w(qi)e"iqi  Rj-3  l  g(qi*q2»q3>UJ(<l2>  l  e-i<13  **3  (13) 

N  *  *  k*i 


k*j 


Expressions  similar  to  (12)  and  (13)  can  be  derived  for  the  case  where  all  ions 
are  not  the  same  by  using,  in  place  of  (10), 


-+  -y 

<k+q|  oj|k>  ~  StqJuj^q)  +  —  \  Ao^(q )e“iq  Ri 


(14) 


where  is  the  host  ion  form  factor  and  is  change  in  form  factor  caused  by 
an  impurity  at  site  R^. 

The  first  term  in  (13)  is  a  local  approximation  to  the  pairwise 

potential  obtained  from  second  order  theory.  The  second  term  is  a  third  order 

correction  which  is  also  pairwise,  and  the  third  term  is  a  sum  over  all  three- 

ion  interactions  involving  sites  i  and  j..  It  is  the  presence  of  this  last 

term  that  makes  V  dependent  on  the  arrangement  of  neighboring  ions.  That  is, 

because  of  the  three-ion  terra,  the  effective  potential  is  not  only  a  function  of 

the  distance  R^j  or  even  the  vector  displacement  Rij,  it  depends  also  on  the 
+  + 

positions  R^  and  Rj  relative  to  neighboring  ions.  As  we  shall  se  in  the  next 
section,  this  has  important  implications  regarding  interatomic  forces  in  the 


vicinity  of  defects 


B.  Applications  to  Defect  Calculations 

In  a  perfect  crystal  the  structure  factor  vanishes  unless  q  equals  a 
recpirocal  lattice  vector  K,  and  all  form  factors  in  (12)  are  evaluated  at 
q  =  K,  where  they  are  quite  weak  for  most  nontransition  metals.  This  is  why  a 
second  order  treatment  is  often  adequate  for  the  calculation  of  perfect  crystal 
properties,  though  significant  third  order  effects  have  been  noted  in 
calculations  of  phonon  spectra  (Ref.  5)  and  cohesive  energies  (Ref.  8)  if 
polyvalent  metals. 

The  structure  factor  is  usually  nonvanishing  at  all  q  in  imperfect 
crystals  and  at  most  q  in  partially  ordered  systems  such  as  liquid  metals,  and 
this  brings  larger,  long  wavelength  values  of  u>(q)  into  play.  Because  the  per¬ 
turbation  is  now  stronger,  one  expects  the  series  (2)  to  converge  more  slowly, 
meaning  that  higher-order,  multi- ion  terms  in  (9)  become  more  important.  In 
liquid  sodium  and  potassium,  for  example,  Hasegawa  (Ref.  9)  found  that  three- 
ion  interactions  result  in  strong,  short  range  attractive  forces,  while  our  own 
calculations  (to  be  described  in  Section  C)  indicate  that  similar  effects  occur 
in  aluminum.  In  their  calculations  of  vacancy  formation  energies  for  sodium, 
magnesium  and  aluminum.  So  and  Woo  (Ref.  10)  found  that,  although  third  order 
terms  are  small,  indicating  reasonably  good  convergence  even  in  this  strong  per¬ 
turbation  case,  the  third  order  contribution  to  the  formation  energy  is  impor¬ 
tant  because  other  terms  nearly  cancel.  For  aluminum,  in  particular,  third 
order  terms  are  essential  to  obtaining  a  positive  formation  energy,  Though 
agreement  with  experiment  is  poor,  the  calculations  still  serve  as  convincing 
evidence  that  third  order,  multi-ion  effects  are  important  in  defect  calcula¬ 
tions. 


More  such  evidence  is  presented  by  Jacucci,  et  al  (Ref.  11)  who  also 
calculated  the  vacancy  formation  energy  as  well  as  the  phonon  spectrum,  liquid 
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structure  factor  and  vacancy  migration  energy  for  aluminum.  Their  calculations 
were  molecular  dynamics  simulations  based  on  a  pairwise  interatomic  potential 
derived  by  Rasolt  and  Taylor  (RT)  (Ref.  12).  The  RT  potential  was  obtained  from 
a  density  functional  calculation  for  an  aluminum  ion  in  jellium,  and  therefore 
contains,  in  effect,  all  of  the  pairwise  corrections  appearing  in  the  perturba¬ 
tion  expansion  of  V^j  in  (9);  it  specifically  excludes  multi- ion  interactions 
involving  three  or  more  ion  sites.  Jacucci,  et  al,  show  that  this  potential 
leads  to  good  agreement  with  all  properties  except  the  vacancy  formation  energy, 
and  argue  that  the  failure  of  this  calculation  is  caused  by  the  exclusion  of 
three- ion  interactions  involving  the  vacancy.  They  conclude  that  one  cannot  use 
a  simple  pairwise  potential  derived  from  bulk  crystal  data,  such  as  phonon  spec¬ 
tra,  to  predict  certain  defect  properties,  because  localized,  multi-ion  intera¬ 
ctions  involving  the  defect  are  not  accounted  for  in  bulk  crystal  data.  Con¬ 
versely,  if  one  uses  experimental  point  defect  data  to  derive  a  pairwise  poten¬ 
tial,  this  potential  should  not  be  used  to  calculate  bulk  crystal  properties 
because  it  does  include  the  effects  of  localized  multi-ion  interactions  at  the 
defect  site. 

This,  we  think,  is  an  extremely  important  point.  It  confirms  an 
assertion  we  made  in  our  proposals  to  ARO,  and  has  been  the  underlying  motiva¬ 
tion  for  our  work  on  multi-ion  interactions.  It  is  important  in  a  basic  sense 
because  the  demonstrated  importance  of  multi-ion  interactions  involving  defects 
refutes  the  idea  that  there  exists  a  universal,  central  pairwise  potential.  It 
is  also  important  technologically,  because  the  central  potential  assumption 
underlies  almost  all  studies  of  defect  properties  by  the  molecular  dynamics 
method,  the  lattice  statics  method  and  other  techniques  for  predicting  the 
structures  and  interacitons  of  crystalline  defects.  If  one  cannot  use  the  same 
potential  in  a  defect-free  region  as  one  uses  near  a  defect,  then  there  is 


serious  reason  to  question  the  validity  of  such  defect  studies. 

Our  research  has  been  concerned  with  simple  defects  in  aluminum.  We 
have  chosen  aluminum  because  it  is  amenable  to  a  pseudopotential/perturbation 
treatment,  its  perfect  crystal  characteristics  are  well  understood,  and  because 
experimental  data  on  its  defect  properties  are  readily  available.  Also,  because 
aluminum  is  a  polyvalent  metal,  we  expect  third  order  effects  to  be  more  promi¬ 
nent  than  they  would  be,  for  example,  in  the  alkali  metals  (Ref.  8,10). 

We  have  been  concerned  mostly  with  third  order  effects,  though  some  of 
our  earlier  work  dealt  with  second  order  calculations  (Ref.  3)  and  with  a 
t-matrix  version  of  the  multi- ion  expansion  (Ref.  4)  for  use  with  strong 
scatterers.  We  are  interested  in  determining  the  circumstances  under  which  a 
third  order  treatment  is  sufficient,  and  in  exploring  some  of  the  properties  of 
third  order  forces.  Because  a  nonlocal  third  order  calculation  is  not  feasible, 
we  have  also  studied  the  use  of  various  local  approximations  to  model  potential 
form  factors  for  use  in  third  order  terms.  More  recently,  we  completed  a  third- 
order  .attice  statics  calculation  of  the  lattice  displacement  field  around  point 
defects  following  the  theory  of  Solt  and  Werner  (Ref.  14).  Tftis  work  is  summa- 
arized  in  the  next  section. 

C.  Summary  of  Progress 
1 .  Early  Work 

Our  first  defect  studies  were  second  order  calculations  of  stack¬ 
ing  fault  energies  for  various  nontransition  metals  (Ref.  13,  Appendix  A).  This 
is  a  situation  where  a  second  order  treatment  should  suffice  because  structure 
factors  for  stacking  faults  are  such  that  the  large  values  of  atomic  form  fac¬ 
tors  at  small  q  do  not  enter  the  calculation  (Ref.  1,2).  Agreement  with  experi¬ 
ment  was  reasonable  considering  the  uncertainties  in  pseudopotentials,  screening 
fucntions  and  experimental  data  that  existed  at  that  time. 
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In  considering  applications  to  other  defects  and  to  transition 
metals/  however,  it  was  clear  that  the  strong  perturbations  encountered  in  such 
cases  required  a  higher  order  treatment.  We  therefore  developed  a  t-matrix 
generalization  of  Harrison's  multi-ion  theory  (Ref.  3),  and  demonstrated  that 
three-ion  interaction  energies  are  indeed  quite  large  for  transition  metals 
(Ref.  4,  Appendix  B) . 

An  application  to  stacking  faults  in  transition  metals  followed 
(Ref.  15,  Appendix  C) .  We  used  an  approximate  version  of  the  third  order  t- 
matrix  theory  and  obtained  reasonable  agreement  with  experimentally  determined 
stacking  fault  energies.  Predictions  of  the  relative  stabilities  of  close- 
packed  structure  were  obtained  from  these  calculations  and  agreement  with  ob¬ 
served  stable  structurs  was  obtained  in  almost  all  cases.  A  somewhat  surprising 
result  was  that  resonant  scattering  of  d-electrons  made  three-ion  interactions 
dominant  in  transition  metals,  these  energies  being  about  ten  times  as  large  as 
pairwise  contributions  to  the  stacking  fault  energy. 

We  then  returned  to  the  study  of  defects  in  nontransition  metals 
and  attempted  to  formulate  a  third  order  treatment  of  defect  scattering  in  terms 
of  the  electron  Green  function  for  the  perfect  crystal.  We  used  a  first  order, 
degenerate  perturbation  theory  calculation  of  the  Green  function,  because  we 
were  concerned  about  the  influence  of  discontinuities  at  Brillioun  zone  bound¬ 
aries  on  defect  screening  charge  distributions.  Although  the  results  are  quite 
complicated,  we  were  successful  in  developing  a  spherical  harmonics  expansion  of 
this  Green  function  for  use  in  a  partial  wave  treatment  of  point  defect  scat¬ 
tering.  Unfortunately,  the  large  number  of  terms  (as  many  as  50)  needed  in  the 
expansion  made  such  a  calculation  impractical,  and  the  idea  of  performing  a  full 
calculation  was  abandoned.  Instead  we  applied  the  Green  function  to  an  ideal¬ 
ized  delta  function  scatterer  and  calculated  the  resulting  screening  charge 


density,  and  host-ion  defect  interaction  energy.  The  results,  which  were  not 
published,  indicated  that  host  crystal  scattering  is  important  and  that  it  leads 
to  strong  anisotropy  in  both  the  screening  charge  density  and  the  host-ion  de¬ 
fect  interaction  energy. 

To  permit  the  use  of  more  realistic  defect  potentials,  we  then 
turned  to  a  third-order,  nondegenerate  perturbation  treatment,  as  given  by  Lloyd 
and  Sholl  (Ref.  6),  and  Brovman  and  Kagan  (Ref.  7,16).  We  chose  to  concentrate 
on  second-order  calculations  of  the  point  defect  screening  charge  density  be¬ 
cause  such  calculations  are  formally  similar  to  third  order  calculations  of  the 
effective  interatomic  potential,  and  both  calculations  involve  the  same  physical 
assumption,  namely,  that  second  order  screening  is  adequate. 

Another  reason  for  concentrating  on  screening  charge  densities  for 
point  defects  is  that,  according  to  Kohn  and  Vosko  (Ref.  17)  and  Blandin  and 
Friedel  (Ref.  18),  the  electric  field  gradient  (EFG)  at  a  host  ion  site  a  dis¬ 
tance  R  from  an  impurity  ion  is  given  by 

q(R)  =  -  a  6ps(R)  (15) 

where  e  is  the  electronic  charge,  5ps(R)  is  the  change  in  the  smooth  part  of  the 
valence  wave  function  (that  part  determined  by  pseudopotential/perturbation 
theory)  and  a  is  a  dimensionless  constant  that  corrects  for  the  effect  of  the 
oscillatory  part  of  the  valence  electron  density  inside  the  host  ion  core.  The 
EFG,  which  can  be  determined  from  nuclear  magnetic  resonance  experiments,  is 
therefore  proportinal  to  6ps,  which  is  calculated  quantity,  thus  providing  a 
basis  for  verifying  screening  charge  calculations.  The  fact  that  our  prelimi¬ 
nary  calculations  indicated  a  strong  host  crystal  scattering  effect  on  6ps  made 
such  comparisons  all  the  more  interesting,  because  experimental  data  on  the 
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EFG's  at  second  neighbor  sites  in  aluminum  alloys  are  unusually  small,  a  fact 
that  had  been  attributed  to  some  unknown  influence  of  the  host  crystal  (Ref.  19). 
Thus,  the  idea  was  that  a  second  order  calculation  of  6ps  that  includes  host 
crystal  scattering  might  serve  to  explain  the  second  neighbor  EFG  anomaly  in 
aluminum. 

Preliminary  second  order  calculations  were  performed  for  five  im¬ 
purities  using  damped  Ashcroft  potentials  (Ref.  20)  with  parameters  adjusted  to 
obtain  rough  agreement  with  the  energy-wavenumber  data  of  Appapillai  and  Williams 
(Ref.  21).  The  results  once  again  showed  strong  anisotropy  in  6ps  with  a  tenden¬ 
cy  for  first  and  second  order  terms  to  cancel  at  the  second  site.  While  this 
aspect  of  the  calculation  was  encouraging,  overall  agreement  with  experiment  was 
quite  poor. 

This  led  us  to  consider  several  possible  improvements  in  the  calcu¬ 
lations.  The  first,  and  most  seriuos  question  raised  by  the  results  concerned 
the  adequacy  of  a  second-order  treatment.  The  fact  that  second  order  terms  were 
large  seemed  to  say  that  third  and  possibly  even  higher  order  corrections  are 
also  important.  On  the  other  hand,  the  pseudopotentials  used  were  rather  crudely 
constructed,  and  we  could  not  be  sure  that  they  were  realistic.  Also,  because 
were  were  concentrating  on  anisotropic  effects  involving  scattering  by  the  defect 
and  host  crystal,  we  did  not  include  the  isotropic  second  order  term  corre¬ 
sponding  to  double  scattering  at  the  defect  site.  Another  possible  source  of 
error  was  that  lattice  distortion  effects  were  completely  ignored.  Finally,  from 
the  derivation  of  (15)  it  is  clear  that  if  one  carries  the  calculation  of  6ps  to 
second  order,  then  the  formula  for  the  core  enhancement  factor  should  be  modified 
accoridngly,  and  this  was  not  done  in  our  calculations.  Attempts  to  incorporate 
these  improvements  in  the  calculations  are  described  below. 
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2.  Nonlinear  Screening  of  Point  Defects 
a.  Je Ilium  Model 

As  a  first  test  of  a  second  order  theory,  screening  charge 
densities  were  calculated  for  the  proton,  an  isolated  aluminum  ion,  and  for 
substitutional  magnesium  and  zinc  in  a  neutral,  uniform  medium  (jellium)  with 
electron  density  equal  to  that  in  the  aluminum  crystal  (Appendix  F) .  To 
calculate  the  density  to  second  order  we  use  the  following  expressions,  which 
are  derived  in  Appendix  F: 


•*<«>  -  £  i£SSHb*«>  * 


(16) 


where  u(q) ,  the  form  factor  self  consistent  to  second  order,  satisfies, 

u)(g)  =  ojjutq)  +  i'^qft  N  \  9<q»P»r>  w(p)  u>(r)  (17) 

P 
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where  wj^(q)  =*  w^on(q)/e(q)  is  the  linear  screening  form  factor  and  r  =  -(q+p) . 

An  approximate  solution  for  w(q) ,  which  is  accurate  through  third  order  in  w^(q) 
(Appendix  E),  is  obtained  by  substituting  u)£(q)  for  u(q)  inside  the  sum;  succes¬ 
sively  better  approximations  can  then  be  generated  by  repeated  substitutions. 
Self  consistent  numerical  results  obtained  by  this  iterative  method,  by  the 
approximate  solution  for  u(q) ,  and  by  linear  screening  theory  are  presented  in 
Appendix  F. 

For  a  substitutional  defect,  the  perturbation  is  taken  as  the 
difference  between  the  impurity  and  aluminum  potentials,  and  is  therefore  much 
weaker  than  an  isolated  ion  potential.  The  results  show  that  for  substitutional 
defects  second  order  effects  are  small,  and  that  the  screening  charge  density 
obtained  by  iteration  to  full  self  consistency  differs  very  little  from  that 
obtained  with  the  approximate  solution  for  u>(q) .  As  one  might  expect,  the 
stronger  perturbation  associated  with  the  aluminum  ion  potential  (or,  with  a 
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change  of  sign,  the  aluminum  vacancy  potential)  produces  somewhat  greater  second 
order  effects  and  Is  more  dependent  on  approximations  to  self  consistency, 
while  the  strongest  scatterer,  the  proton,  shows  very  large  second  order  effects 
and  still  greater  dependence  on  self  consistency. 

Comparison  of  the  proton  calculations  with  exact  density 
functional  results  (Ref.  22)  shows  that  adding  second  order  terms  gives  much 
better  results  but  does  not  fully  account  for  nonlinear  screening  effects.  For 
the  aluminum  ion  there  is  some  disagreement  between  published  density  calcula¬ 
tions  (Ref.  23,24)  and  comparisons  of  second  order  calculations  with  these 
results  is  therefore  inconclusive.  However,  it  is  encouraging  to  note  that  the 
disagreement  between  second  order  and  density  functional  results  is  no  greater 
than  the  disagreement  among  the  density  functional  calculations  themselves. 

For  substitutional  impurities,  the  fact  that  nonlinear  screening  corrections 
are  small  suggests  that  a  second  order  theory  of  the  screening  charge  density 
is  adequate  for  such  defects. 

b.  Host  Crystal  Effects 

When  host  ion  scattering  is  added  to  the  jellium  model,  the 
Fourier  transform  of  the  screening  charge  density  is  still  given  by  (16),  but 
with  anisotropic  w(q)  and  6ps(q) .  We  have  not  attempted  a  fully  self  consistent 
treatment  in  this  case  because  of  the  computational  difficulty  it  presents. 
Instead  we  calculated  the  self  consistent  spherical  average  of  u>(q),  which 
satisfies 

u(q)  -  w°(q)  +  t^q?*  ^  I  <*£  7^  Jg(q*K*p) w(p)dQp  (18) 
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where  is  the  host  crystal  form  factor  at  q  ■  K. 
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Hie  screening  charge  density  Is  then 
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Typical  results  for  substitutional  magnesium  are  presented  in 
Appendix  G,  in  this  case  for  average  local  potentials  based  on  damped  Heine- 
Abarenkov  model  potentials  for  magnesium  and  aluminum,  these  data  show,  once 
again,  that  anisotropic  contributions  to  the  screening  charge  density  are 
important,  though  they  are  not  as  pronounced  as  our  earlier  calculation  with 
delta  function  potentials  would  indicate.  As  it  turns  out,  however,  second 
order  corrections  to  the  EFG  are  much  smaller  than  one  would  expect,  and  for 
this  reason,  calculations  of  anisotropic  screening  were  not  carried  beyond  this 
point.  The  role  of  higher  order  terms  in  the  EFG  formula  is  discussed  below. 

3.  Electric  Field  Gradient  Calculations 
a.  Second  order  Effects 

If  we  let  6p(r)  be  the  total  change  in  valence  electron 
density  (including  host  ion  oscillations)  then  the  radial  component  of  the 
valence  effect  EFG  (the  EFG  in  the  absence  of  lattice  distortion)  is  given  by 


(Ref.  25) 
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where  P2(r*R)  is  a  Lengendre  polynominal,  and  Az  is  the  impurity  valence  minus 
the  host  valence. 
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Following  Kohn  and  Vosko  (Ref.  17)  we  might  now  assume  that 
+  * 

6p(q)  can  be  replaced  by  a 6pg(q)  where  a  is  the  core  enhancement  factor.  How¬ 
ever,  as  we  show  in  Ref.  26  (Appendex  E),  if  we  include  anisotropic  second  order 
terms  in  6ps(q) ,  then  such  terms  should  also  be  included  in  the  calculation  of 
a.  Hiis,  it  turns  out,  leads  to  anisotropic  corrections  to  a  of  the  order  of 
20%. 

We  might  also  assume  that  R  is  large  and  that  the  only  impor- 

■> 

tant  q  directions  are  therefore  those  that  are  nearly  parallel  to  R,  in  which 
case  p2<q.R)~1,  and  the  last  term,  as  well  as  the  second  term,  becomes  propor¬ 
tional  to  6ps.  However,  in  our  calculations  we  chose  to  avoid  this  particular 
asymptotic  approximation  because  numerical  evaluation  of  the  q  integral  is 
required  in  any  case. 

Typical  results  are  reported  in  Appendix  G.  Hie  two  principal 
points  to  be  noted  here  are  ( 1 )  the  linear  screening  part  of  the  last  term  in 
(21)  is  the  dominant  contribution  at  R  of  interest,  and  (2)  this  term  is  not 
proportional  to  6p„(R),  contrary  to  the  Kohn-Vosko  formula. 

This  totally  unexpected  result  tells  us  that  the  long  standing 
disagreement  between  calculated  and  measured  EFG's  is  caused  not  so  much  by 
inaccuracies  in  the  calculation  of  6ps  or  a,  as  by  the  failure  of  the  Kohn-Vosko 

A  A 

formula  itself.  In  particular,  the  asymptotic  argument  that  leads  to  P2(q*R)~1 
produces  serious  errors  at  all  R  of  interest. 

This  being  the  case,  we  were  led  to  question  the  validity  of 
the  core  enhancement  approximation,  6p(R)  ~  a6ps(R),  because  this  approximation 
is  also  based  on  a  large  R  assumption.  For  this  reason  we  decided  to  abandon 
the  Kohn-Vosko  formula,  return  to  the  original  EFG  expression  (21),  and  recal¬ 
culate  the  EFG  without  invoking  an  asymptotic  approximation.  These  calculations, 
which  are  reported  in  Appendix  G,  were  done  in  the  linear  screening  approximation 
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because  Table  I  of  Appendix  G  shows  that  second-order  effects  are  very  small, 
b.  A  Correction  to  the  Kohn-Vosko  Formula 

As  reported  in  Appendix  G,  we  used  the  core  state  projection 
operator  technique  (Ref.  1)  to  calculate  the  impurity-induced  change  in  the 
oscillatory  part  of  the  valence  wave  function  in  host  ion  cores.  This  result 
was  then  used  in  (21)  to  calculate  the  core  enhancement  part  of  the  EFG.  Not 
surprisingly,  the  numerical  results  show  that  core  enhancement  effects  cannot  be 
accounted  for  by  a  simple  core  enhancement  factor  that  is  a  property  of  the  host 
crystal  alone. 

On  the  other  hand,  we  also  show  in  Appendix  G  that  the  XOhn- 
Vosko  EFG  formula  results  if  we  make  the  approximation 

j2(qR)  ~  j0(qR)  (22) 

which  holds  for  large  qR,  and  also  assume  that  certain  core  state  matrix  ele¬ 
ments  that  appear  inside  the  q  integral  can  be  approximated  by  their  values  at 
the  singular  point  q  =  2fcF.  This,  as  we  said,  leads  to  the  Kohn-Vosko  approxi¬ 
mation  and  gives  a  =  6.8,  in  reasonable  agreement  with  more  elaborate  calcula¬ 
tions  (Ref.  27). 

The  reason  this  approximation  fails  is  that  (22)  is  not  valid 
at  small  q  and,  when  Az*0,  the  small  q  part  of  the  EFG  integrand  is  quite 
important.  To  demonstrate  this  we  calculated  a  correction  term  to  account  for 
the  error  in  (22)  and  found  that  when  Az*0,  this  new  term  dominates  at  all  R. 
Finally,  returning  to  the  large  R  assumption,  we  found  the  following  corrected 
version  of  (15): 


where  the  first  term  is  the  correction  for  small  q  (long  wavelength)  effects  and 
F(q)  is  a  function  of  core  state  matrix  elements.  Data  presented  in  Appendix  G 


show  that  the  new  asymptotic  formula  gives  EFG  predictions  in  good  agreement 
with  a  full  nvxnerical  treatment  even  at  nearest  neighbor  sites. 

4.  Lattice  Distortion  Calculations 

In  Section  B  we  cited  reasons  why  we  believe  a  point  defect/ 
lattice  distortion  calculation  should  be  carried  to  third  order.  To  repeat,  we 
agree  with  Jacucci,  et  al  (Ref.  11),  that  their  calculations  demonstrate  that 
three-ion  interactions  are  important  in  point  defect  studies  such  as  theirs,  and 
that  some  such  interactions  are  not  accounted  for  when  one  uses  a  semiempirical 
interatomic  potential  derived  from  bulk  crystal  data.  Furthermore,  according  to 
Solt  (Ref.  28),  the  asymptotic  expression  for  the  displacement  field,  as  calcu¬ 
lated  by  lattice  statics  (Refs.  29,30)  is  incorrect  unless  one  includes  third 
order  terms. 

Therefore,  encouraged  by  the  results  of  our  second  order  calcu¬ 
lations  of  defect  screening  charge  densities  (Appendix  F) ,  we  undertook  a  third 
order,  lattice  statics  calculation  of  the  ion  displacement  field  caused  by  a 
vacancy  in  aluminum.  The  expression  we  used  for  the  third  order  force  was  that 
given  by  Solt  &  Werner  (Ref.  14),  while  the  perfect  crystal  dynamical  matrix  was 
calculated  to  second  order  using  the  empirical  model  potential  of  Sen,  et  al 
(Ref.  31). 

The  results  are  reported  in  Appendix  H.  They  show  that,  in  second 
order,  our  calculated  displacements  are  significantly  smaller  than  those  ob¬ 
tained  by  Singhal  in  an  earlier  calculation  (Ref.  32).  This  is  not  particularly 
surprising  in  view  of  the  fact  that  Singhal ’#  model  potential  differs  from  that 
of  Sen,  et  al.  However,  when  third  order  forces  were  added,  the  displacements 
became  still  smaller,  in  some  cases  even  changing  sign,  which  indicates  that 


second  and  third  order  forces  are  of  the  same  order  of  magnitude 
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This  raises  some  serious  questions  concerning  the  use  of  third  order  perturba¬ 
tion  theory  in  the  calculation  of  defect  properties.  In  particular,  if  we  are 
to  believe  that  third  order  forces  are  as  strong  as  our  calculation  shows,  then 
we  might  expect  that  higher  order  multi-ion  interactions  are  also  strong  and 
cannot  be  ignored.  On  the  other  hand,  perhaps  our  use  of  the  local  approxima¬ 
tion  tends  to  exaggerate  third  order  forces  thus  leading  to  unreasonably  large 
third  order  displacements.  In  either  case,  our  results  tend  to  cast  doubt  on 
the  validity  of  a  local,  third  order  theory  of  interionic  forces  involving  de¬ 
fects. 

D.  Conclusions 

For  the  most  part,  our  research  has  focused  on  the  calculation  of 
screening  charge  densities  around  defects  in  nontransition  metals.  This  is 
because  changes  in  the  screening  charge  distribution  induced  by  defects  are  of 
fundamental  importance  in  determinig  interionic  forces  among  neighboring  ions. 

Our  calculations  indicate  that  a  pseudopotential/perturbation  calcu¬ 
lation  carried  to  second  order  is  adequate  for  the  calculation  of  charge  densi¬ 
ties  around  weak  scatterers,  such  as  substitutional  impurities  with  small  valence 
difference.  In  fact,  in  such  cases  we  find  that  a  linear  screening  (first 
order)  treatment  should  suffice  for  most  purposes.  Even  for  strong  scatterers 
like  the  vacancy  in  aluminum,  the  theory  is  in  reasonable  agreement  with  more 
accurate  density  functional  calculations. 

Comparisons  with  experimental  data  on  electric  field  gradients  (EFG) , 
around  point  defects,  which  are  determined  in  part  by  the  charge  density  distri¬ 
bution,  are  inconclusive.  This  is  because  the  calculations  are,  in  a  sense, 
incomplete  in  that  an  accurate  treatment  of  the  lattice  distortion  contribution 
to  the  EFG  is  still  lacking.  Still,  from  the  calculations  themselves  we  can 
draw  two  important  conclusions  regarding  EFG  theory.  These  are  the  following: 

( 1 )  a  first  order  calculation  of  the  screening  charge  density  is  sufficient  for 


EFG  predictions,  and  (2)  the  approximate  Kohn-Vosko  formula,  which  relates  the 
EFG  to  the  screening  charge  density,  misses  an  important  contribution  to  the  EFG 
in  cases  where  the  host  and  impurity  valences  differ.  An  approximate  expression 
for  this  missing  contribution  was  derived  and  is  reported  in  Appendix  G. 

We  conclude,  therefore,  that  second  order  perturbation  theory  gives  a 
good  account  of  defect  screening  in  all  but  a  few  cases  where  the  defect  pseudo¬ 
potential  is  extremely  strong.  In  fact,  from  the  opinions  expressed  in  the 
literature  on  defect  screening,  we  would  say  that  the  theory  is  much  better  than 
expected. 

The  theory  of  interionic  forces  is,  however,  another  matter.  Our 
calculation  of  the  effect  of  three-ion  forces  on  lattice  distortion  around  a 
vacancy  produced  third  order  displacements  that  are,  in  our  opinion,  much  too 
large.  It  is  unclear  at  present  whether  this  was  caused  by  the  omission  of 
higher  order  forces  or  by  the  local  potential  used  in  our  third  order  term. 
Because  both  the  calculation  of  higher  order  terms  and  the  calculation  of  non¬ 
local  effects  in  third  order  are  extremely  difficult,  it  seems  unlikely  that 
this  question  can  be  easily  answered. 

We  are  left,  then,  with  a  strong  indication  that  the  path  we  have  been 
pursuing,  i.e.,  a  low  order,  local  peturbation  theory  of  interionic  forces,  no 
longer  appears  promising.  This  leaves  only  density  functional  methods,  which 
are  so  far  restricted  to  spherically  symmetric  systems,  and  self-consistent 
cluster  tyupe  calculations  as  candidate  methods  for  calculating  interionic 
forces.  Unfortunately,  neither  of  these  approaches  is  well-suited  to  a  force 
calculation  where  defects  are  involved,  and  it  is  therefore  difficult  to  recom¬ 
mend  a  direction  for  further  research.  At  present,  our  personal  choice  would  be 
to  use  density  functional  theory  to  derive  accurate  pairwise  potentials  and  then 
look  for  approximate  ways  to  correct  these  potentials  for  multi-ion  interactions 
associated  with  the  presence  of  neighboring  ions. 
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Calculations  of  stacking-fault  energies  in  aluminum,  magnesium,  beryllium,  copper,  silver,  and  gold  are 
reported.  For  the  polyvalent  metals  it  is  shown,  by  means  of  comparisons  of  numerical  results  based  on 
several  different  energy-wave-number  characteristics,  that  exchange  and  correlation  corrections  are  not 
as  significant  in  the  present  context  as  they  are  in  most  other  defect  calculations.  Nonlocal  effects  are 
quite  important  but  can  be  approximately  accounted  for  by  empirical  adjustments  of  the  form  factors 
based  on  simpler  local  models.  Although  stabilities  against  faulting  and  the  relative  magnitudes  of 
stacking-fault  energies  are  correctly  predicted,  quantitative  agreement  with  experiment  is  not  obtained. 
Possible  reasons  for  this  discrepancy  are  discussed.  Results  for  the  noble  metals  indicate  a  severe 
sensitivity  to  the  detailed  curvature  of  the  energy-wave-number  characteristic  in  the  vicinity  of  the 
smallest  reciprocal-lattice  vector.  Failure  to  obtain  agreement  with  experiment  in  these  cases  may 
therefore  be  due  to  minor  inaccuracies  in  the  energy-wave-number  characteristics  for  noble  metals. 


I.  INTRODUCTION 

Applications  of  the  pseudopotential  theory  of 
metals  to  analyses  of  the  energetics  of  most  crys¬ 
talline  defects  are  made  difficult  by  a  number  of 
factors.  Among  these  are  the  questionable  validity 
of  a  first-order  perturbation  treatment,  1  the  sen¬ 
sitivity  of  interionic  potentials  to  exchange  and 
correlation  corrections,  2,3  and  the  need  to  proper¬ 
ly  account  for  the  nonlocal  nature  of  the  pseudo¬ 
potential  itself.4  However,  the  calculation  of 
stacking-fault  energies  in  close-packed  crystals 
is  one  class  of  problems  that  apparently  does  not 
suffer  from  the  first  mentioned  difficulty,  owing 
essentially,  to  the  fact  that  close  packing  is  pre¬ 
served  across  the  fault  plane. 1  There  is  also 
reason  to  expect  that  stacking-fault-energy  pre¬ 
dictions  are  less  sensitive  than  other  defect  cal¬ 
culations  to  uncertainties  in  exchange  and  corre¬ 
lation  effects.  This  is  because  effective  inter¬ 
ionic  potentials,  which  are  always  used  at  least 
implicitly  in  defect-energy  calculations,  are  most 
sensitive  to  such  many-body  corrections  at  small 
interionic  distances  (or  the  order  of  the  neares' 
neighbor  distance), 3  while  the  smallest  distance 
involved  in  a  stacking-fault-energy  calculation  is 
twice  the  distance  between  close-packed  planes.  3 
The  effects  of  nonlocality,  however,  are  more 
difficult  to  estimate  because  such  effects  persist, 
in  the  form  of  a  reduction  in  magnitude  of  the  os¬ 
cillations  in  the  interatomic  potential,  at  large 
interatomic  distances.4  On  the  other  hand,  fully 
nonlocal  treatments  do  exist,  2,8-8  so  that  it  is 
now  possible  to  avoid  this  difficulty  altogether. 
Thus,  if  our  suggestion  that  the  stacking-fault  en¬ 
ergy  is  insensitive  to  many-body  corrections  is 
valid,  it  should  be  possible  to  proceed  with  rea¬ 
sonable  confidence  in  the  prediction  of  stacking- 
fault  energies  for  at  least  a  few  close-packed 
metals. 


In  the  present  paper  we  examine  the  effects  of 
both  many-body  corrections  and  nonlocality  by 
means  of  comparisons  of  calculated  stacking-fault 
energies.  CXir  principal  purpose  in  undertaking 
this  work  was  to  test  the  sensitivity  of  the  result 
to  such  effects  and  thus  to  investigate  the  useful¬ 
ness  of  the  more  rigorous  models,  and  simple 
local  potentials  as  well,  for  predicting  trends  in 
stacking-fault  energies.  A  secondary  objective 
was  to  explore  possible  advantages  of  the  calcu¬ 
lations!  method  of  Blandin,  Friedel,  and  Saada, 5 
in  which  the  stacking-fault  is  expressed  as  a  real- 
space  sum,  as  opposed  to  the  usual  reciprocal- 
space  formulation,  8  involving  the  effective  interac¬ 
tion  between  close-packed  planes. 

A  review  of  the  method  of  calculation  and  a  de¬ 
scription  of  the  scope  of  the  work  are  presented 
in  Secs.  II  and  in.  In  Sec.  IV  we  present  numeri¬ 
cal  results  for  the  interplanar  potential  in  alumi¬ 
num  and  gold,  these  results  being  typical  of  poly¬ 
valent  and  noble  metals,  respectively.  Topics 
discussed  here  include  the  effects  of  nonlocality 
and  many-body  corrections  in  polyvalent  metals, 
the  asymptotic  approximations  of  Blandin  ei  at. ,  5 
the  sensitivity  of  the  gold  potential  to  minor  fea¬ 
tures  of  the  energy-wave-number  characteristic, 
and  implications  of  this  sensitivity  as  regards  the 
relative  stabilities  of  the  two  close-packed  struc¬ 
tures  for  noble  metals.  Section  V  contains  the  re¬ 
sults  of  stacking-fault-energy  calculations,  com¬ 
parisons  with  experiment,  and  an  analysis  based 
on  the  discussion  in  Secs.  Ill  and  TV. 

II.  METHOD  OF  CALCULATION 

The  stacking- fault  energy,  which  will  be  de¬ 
noted  by  y,  Is  defined  as  the  difference  In  energy 
per  unit  fault  area  between  the  faulted  and  perfect 
crystals.  In  the  usual  pseudopotential- perturbation 
approximation,  the  expression  for  y  is  a  sum  over 
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FIG.  1.  Normalized  cue rgy — wave-nutnbe r  characteris¬ 
tics  for  aluminum.  Stacking-fault  energy  calculations 
involve  only  those  data  for  which  q  is  greater  than  the 
smallest  reciprocal-lattice  vector,  indicated  here  by  the 
line  labeled  y.  Sources  are  as  follows:  A,  Animalu 
(Ref.  121;  Y,  Yamamoto  (Ref,  8);  SP,  Shaw  and  Pynn 
(Kef.  71;  and  EC,  the  empty-core  model  discussed  in  the 
text.  The  energy-wave-number  characteristics  of  Shaw 
(Ref.  6)  and  Shyu  vt  al.  (Ref.  2)  differ  only  slightly  from 
the  Shaw  — Pynn  curve. 


wave-vector  space  of  a  factor  that  depends  only  on 
q,  the  electron  wave  number,  multiplied  by  the 
difference  in  the  absolute  squares  of  the  faulted- 
and  perfect-crystal  structure  factors. 9  Blandin 
ft  al.,1  have  shown  that  this  expression  can  be 
reduced  to  the  simple  form 

S  S(n)Aq>(nh)  ,  (1) 

»»2 

where  h  is  the  distance  between  close-packed 
planes,  ,V(»)  is  a  set  of  integral  weights  corre¬ 
sponding  to  a  particular  fault  configuration,  and 
&<p(nh)  is  an  interplanar  potential  difference  de¬ 
fined  as  follows;  If  the  stacking  sequence  of  close- 
packed  planes  is  denoted  in  the  usual  way  by  a 
sequence  of  symbols  (e.  g. ,  ABABAB. . .  for  the 
hep  structure),  then  A<p(nfi)  is  the  interaction  en¬ 
ergy  of  a  pair  of  planes  in  nonequivalent  (e.g., 
A-B)  positions  minus  the  energy  of  a  pair  of  planes 
in  equivalent  (e.g.,  A-A  or  B-B)  positions,  the 
interplanar  distance  in  both  cases  being  nh.  The 
reason  that  the  sum  in  Eq.  (1)  begins  with  n  =  2  is 
that  in  all  of  the  stacking  sequences  considered 
here,  nearest- neighbor  planes  are  in  nonequiva¬ 
lent  positions  in  both  the  perfect  and  faulted  crys¬ 
tals  and  their  interaction  energy  therefore  cancels 
in  calculating  the  energy  change  due  to  faulting. 

The  function  A<p(nh\  is  given  by 

?  fcos  g.  d-  nh) ,  (2) 


where  z  is  the  effective  valence,  1  kr  is  the  Fermi 
wave  number,  a  is  the  lattice  constant  for  the 
close- packed  plane,  g  is  a  reciprocal  lattice  vec¬ 
tor  for  the  close-packed  plane,  g-  Igl,  and  ±d  is 
tne  parallel  displacement  of  one  plane  relative  to 
a  second  plane  in  a  nonequivalent  position.  1  The 
dimensionless  function  J(g,  nh)  is  the  integral 

J(g,  »*)  =  J^  p  ~ F^i}xt  X~* — co&k rnhr)dt]  , 

(3) 

where  x=g/kr  and  FM{q/kr)  is  the  normalized  en¬ 
ergy-wave-number  characteristic. 4  In  practice 
the  only  significant  terms  in  Eq.  (2)  are  those  for 
the  smallest  nonvanishing  value  of  g.  (By  actual 
numerical  test,  higher-order  terms  contribute  less 
than  1%  to  A q>. )  In  all  subsequent  discussions  we 
will  therefore  use 

AQ  *1 

&q>(nh)= — k^a*’  Jfah)  ,  (4) 

which  is  obtained  by  summing  over  the  six  smallest 
g  vectors.  The  integral  J{nh),  which  is  propor¬ 
tional  to  the  interplanar  potential,  is,  of  course, 
given  by  the  right  side  of  Eq.  (3)  with  x  equal  to  its 
minimum  value.  Actually,  as  the  development  of 
Eq.  (2)  reveals, 1  the  interplanar  potential  is  pro¬ 
portional  to  +J(nh)  for  planes  in  equivalent  positions 
and  ~J(nh)  for  planes  in  nonequivalent  positions. 

From  Eq.  (3)  it  is  evident  that  the  only  values 
of  the  argument  of  F „  that  enter  the  calculation 
are  those  for  which  q/kF  >  x.  In  Figs.  1  and  2  we 
show  the  location  of  x  with  respect  to  the  peaks  and 
valleys  in  the  energy-wave- number  characteristics 


q/kp 


FIG.  2.  Normalized  energy-wave-number  character¬ 
istics  for  gold.  The  line  marked  y  has  the  same  signif¬ 
icance  as  in  Fig.  1.  The  solid  curve  is  based  on  the 
tabulation  of  Morlarty  (Ref,  lu)  ar.d  the  dashed  curve  Is 
an  arbitrary  alteration  of  t^a’  data  as  described  ’ n  the 
text. 
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for  aluminum  and  gold;  these  data  are  typical  of 
polyvalent  and  noble  metals,  respectively.  Since 
X  is  of  the  order  of,  or  greater  than,  the  position 
of  the  first  minimum  in  FN ,  the  large  matrix 
elements  of  the  pseudopotential  that  occur  near 
q/kf=  0  do  not  enter  the  calculation.  Thus,  as 
Heine  and  Weaire1  point  out,  the  first-order  per¬ 
turbation  approximation  on  which  the  formalism 
is  based  is  more  likely  to  be  valid  for  stacking- 
fault  calculations  than  it  is  for  other  defect  studies 
that  do  involve  values  of  q/kr  near  zero. 

Bx<2,  so  that  the  integrand  in  Eq.  (3)  contains 
the  logarithmic  singularity  at  the  point  q=  2kr  (i.e. , 
fj2 »  x*  =  4),  we  can  expect  J(nh)  to  exhibit  a  long- 
range  oscillatory  behavior  similar  to  that  of  inter- 
ionic  potentials.  *  Blandin  et  al. , 5  showed  that,  in 
this  case,  the  asymptotic  form  for  J(nh)  is  (with 
the  usual  assumption  of  a  local  pseudopotential) 

.  (si 

ico  jr  z  n 

where 

v  =  [1  -  (x/2)2] 1/2  (6) 

and  w( 2kr)  is  the  screened-pseudopotential  form 
factor  evaluated  at  q  -  2k  r .  They  also  found  that 
if  this  asymptotic  approximation  is  used  to  calcu¬ 
late  £k<p(nh ),  the  stacking-fault-energy  sum  given 
by  Eq.  (1)  can  be  evaluated  exactly  for  intrinsic, 
extrinsic,  and  twin  faults  in  both  fee  and  hep  crys¬ 
tals.  For  x  >2  the  poles  in  Eq.  (3)  occur  at  imag¬ 
inary  J)  and  the  resulting  asymptotic  form  is  a 
'ecaying  exponential.  In  this  case  the  sura  in  Eq. 
vl)  converges  quite  rapidly. 

In  the  present  case,  for  x  <  2,  we  used  the  asymp¬ 
totic  formulas  of  Blandin  et  al.,  with  numerically 
computed  corrections  to  a q>(nh)  for  small  «.  Thus, 
if  we  let  yA  ar>d  &<pA(nli)  denote  the  asymptotic  ap¬ 
proximations  described  above,  then  the  formula 
used  in  place  of  Eq.  (1)  is 

M 

y=yA+  £  !V(n)[A<p(Mfc)-  A(pA(nh)\  ,  (7) 

n>2 

where  M  is  a  number  large  enough  that  the  differ¬ 
ence  between  A(p(Afl;)  and  A <pA(Mh)  is  small. 

As  it  turned  out,  the  most  serious  difficulty  we 
encountered  in  applying  Eq.  (7)  was  in  the  numeri¬ 
cal  evaluation  of  the  integral  J(nh)  for  large  nh. 
Since  preliminary  tests  showed  that  the  mesh -point 
spaefng  in  tabulated  values  of  F„  is  too  coarse  for 
the  application  of  standard  numerical- integration 
methods,  it  was  necessary  to  interpolate  the  data 
to  obtain  estimates  of  Fx  at  intermediate  values  of 
q/k r .  After  experimenting  with  various  interpola¬ 
tion  and  integration  schemes  we  finally  settled  on 
the  Aitken  interpolation  method10  for  obtaining  val¬ 
ues  of  F M  at  ten  points  within  each  tabulated  inter¬ 
val.  The  integration  over  each  subinterval  was 


then  done  by  fitting  the  quantity  in  brackets  in  Eq. 
(3)  with  a  single  exponential  in  >)  and  evaluating 
the  resulting  integral  exactly. 

In  an  attempt  to  minimize  the  inevitable  inaccu¬ 
racy  of  interpolated  values  near  the  logarithmic 
singularity  at  q  ^  2k  t ,  we  performed  the  interpo¬ 
lation  on  the  function  c(q)F N(q/kr)/  (c(q)  -  1],  where 
€(<?)  is  the  Hartree  dielectric  function.  The  idea 
here  was  that  in  the  local  approximation  with 
Hartree  screening  this  function  is  smooth  and, 
therefore,  more  easily  interpolated  than  Fs  itself. 
Even  so,  when  we  tested  the  method  by  applying  it 
to  such  approximate  functions  we  were  not  able  to 
consistently  reproduce  the  asymptotic  oscillations 
in  Jlfth)  for  n  greater  than  about  4.  This  is  illus¬ 
trated  in  Fig.  3,  where  we  show  the  asymptotic 
form  of  J(nh),  the  curve  obtained  by  evaluating  F„ 
exactly  at  each  mesh  point,  and  the  data  obtained 
with  an  interpolated  function.  For  the  test  shown 
here  we  used  an  exponentially  damped,  empty- 
core- model  potential,  11  with  parameters  chosen 
to  obtain  an  approximate  fit  to  tabulated  F*  data  for 
magnesium.  The  lengths  of  the  vertical  lines 
drawn  at  integral  values  of  n  are  proportional  to 
the  weights  -N{n)  for  an  intrinsic  fault  in  an  hep 
crystal.  [From  Eqs.  (1)  and  (3)  it  follows  that  y 
is  proportional  to  the  sum  of  the  products  -  ,V  [n t 
xj(nk).  ] 

From  the  erratic  behavior  obtained  with  the  in¬ 
terpolated  Fn  at  large  values  of  n,  it  is  evident 
that  when  n  is  large  the  numerical  integration  is 
too  sensitive  to  interpolation  errors  to  be  reliable. 

It  is  also  evident  that  the  interpolation  procedure 
leadB  to  some  error  for  small  n,  as  evidenced  by 
results  for  n  between  2  and  3.  These  errors,  how¬ 
ever,  are  not  as  important  because  the  weights 
used  in  computing  y  [the  N{n)  in  Eq.  (7)]  are  greater 


a  (NUMBER  OF  INTERPLANAA  SPAC'NGS) 


FIG.  3.  Interplanar  potential  function  for  magnesium 
based  on  the  local  empty-core  model  described  in  the 
text.  The  lengths  of  the  vertical  bars  are  proportional 
to  the  weights  involved  in  Eq.  (1)  in  the  calculation  o'  the 
intrinsic  fault  energy  in  an  hco  crystal. 
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for  large  n.  Thus,  if  we  use  the  data  shown  in 
Fig.  3  to  compute  the  energy  of  an  intrinsic  fault 
in  the  hep  structure,  we  obtain  11.  5  erg/cm2  with 
the  computed  F*  and  13.  2  erg/cm2  with  the  inter¬ 
polated  function,  if  the  sum  in  Eq.  (7)  is  termi¬ 
nated  at  Af  =  4.  With  M  =  6,  however,  we  obtain 
11.9  and  7.7  erg/cm2  using  the  computed  and  in¬ 
terpolated  functions,  respectively.  Similar  re¬ 
sults  were  obtained  with  empty-core  models  for 
aluminum  and  beryllium,  although  the  errors  in 
y  for  these  data  were  much  smaller. 

Since  interpolation  difficulties  of  this  type  un¬ 
doubtedly  exist  with  tabulated  nonlocal  F„,  the  sum 
in  Eq.  (7)  was  terminated  at  Af  =  4  in  all  of  the  cal¬ 
culations  reported  here.  Fortunately,  this  caused 
no  serious  problem  in  most  of  the  computations  be¬ 
cause  the  agreement  between  computed  and  asymp¬ 
totic  approximations  to  J  was  generally  quite  good 
at  n  =  4.  With  some  of  the  data  for  beryllium,  how¬ 
ever,  the  difference  between  the  computed  and 
asymptotic  J  was  still  significant  at  n  =4  and  the 
fact  that  we  could  not  handle  larger  interplanar  dis¬ 
tances  prevented  us  from  obtaining  reliable  results. 

111.  SCOPE  OF  THE  CALCULATIONS 

The  remaining  sections  of  this  paper  are  con¬ 
cerned  with  the  results  of  calculations  for  alumi¬ 
num,  magnesium,  beryllium,  copper,  silver,  and 
gold.  In  all  cases  stacking-fault  energies  were 
computed  for  intrinsic,  extrinsic,  and  twin  faults 
in  both  the  fee  and  hep  structures,  thus  providing 
tests  of  the  ability  of  the  method  to  predict  stabili¬ 
ty  against  faulting  and  the  relative  energies  of 
different  fault  configurations. 

The  energy-wave-number  data  used  in  calcula- 


FIG.  4.  Interplanar  potential  functions  for  aluminum 
based  on  the  optimized  model  potential.  The  vertical 
bars  correspond  to  weights  for  an  fee  intrinsic  fault. 


tions  for  the  polyvalent  metals  are  indicated  in 
Table  I.  In  comparing  the  results  of  such  calcu¬ 
lations  it  is  important  to  recall  that  there  is  an 
inherent  arbitrariness  in  the  construction  of  a 
model  potential  and  that  this  arbitrariness  leads, 
in  a  truncated  perturbation  expansion,  to  some 
unknown,  model-sensitive  error.  Thus,  as  far 
as  tests  of  sensitivity  to  many-body  corrections 
are  concerned,  the  most  meaningful  comparisons 
we  can  make  are  those  involving  the  first  three 
sets  of  data  listed  in  Table  I.  AH  three  sets  were 
derived  from  the  same  ionic-model  potentials  and 
differ  only  in  the  treatment  of  electron  exchange 
and  correlation  energies.  Similarly,  because  they 
differ  only  in  the  way  nonlocality  is  handled,  the 


TABLE  I.  Energy-wave-number  characteristics  used  in  calculations  for  polyvalent  metals. 


Data  source 

Model  potential 

Screening 

approximation 

calculation 

Shaw-Pynn 
(Ref.  7) 

Optimized* 

Shaw-Pynn 
(Ref.  7) 

nonlocal 

Shaw 
(Ref.  6) 

Optimized 

Hartree 

nonlocal 

Shyu  et  at. 

(Ref.  2) 

Optimized 

SSTLe 

nonlocal 

Animaiu 
(Ref.  12) 

Heine-Abarenkov* 

Hubbard -Sham 11 

semi  local 

Y  amamoto 
(Ref.  8) 

Heine-Abarenkov 

Hubbard -Sham 

nonlocal 

Empty-core 

model 

Exponentially  damped 
empty  core 

Hartree 

local 

*R.  W.  Shaw,  Jr.  and  W,  A.  Harrison,  Phys.  Rev.  163  ,  604  (1967). 

’’References  13-15. 

CK,  S.  Singwi,  A.  SJdlander,  M.  P.  Tosi,  and  R.  H.  Land,  Phys.  Rev.  B  1 ,  1044  (1970). 
llL.  J.  Sham,  Proc.  R.  Soc.  Lond.  A  283,  33  (1965), 
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next  two  sets  uf  data,  those  based  on  the  Heine- 
Abarenkov- model  potential,  provide  a  basis  for 
investigating  the  importance  of  nonlocal  effects. 

The  primary  reason  for  including  calculations 
based  on  the  empty-core  model  was  to  determine, 
through  comparisons  with  results  based  on  the 
more  refined  models,  to  what  extent  the  simpler 
theory  might  be  useful  for  estimating  trends  in 
stacking- fault  energies. 

The  only  nonlocal  energy-wave-number  charac¬ 
teristics  available  for  the  noble  metals  are  those 
reported  by  Moriarty.  16  Although  comparisons 
with  results  based  on  simpler  local  models  could 
have  been  made,  we  did  not  do  this  because  such 
models  are  obviously  unrealistic  for  the  noble  met¬ 
als,  where  the  energy- wave- number  characteris¬ 
tic  is  dominated  by  the  highly  nonlocal  effects  of 
resonant  scattering.  Our  analysis  for  the  noble 
metals  is  therefore  limited  to  comparisons  with 
experiment  and  a  numerical  test  of  the  sensitivi¬ 
ties  of  the  interplanar  potential  and  stacking-fault 
energies  to  minor  alterations  of  the  energy-wave¬ 
number  characteristic. 

IV  INTERPLANAR  POTENTIALS 

The  interplanar  potential  function  J{nh)  for  alu¬ 
minum,  as  detcrnr.ned  from  the  energy-wave¬ 
number  tabulations  of  Shaw6  (Hartree  approxima¬ 
tion),  Shaw  and  Pynn,  ?  and  Shyu  cl  al. , 2  is  plotted 
in  Fig.  4  along  with  the  asymptotic  approximation 
based  on  the  Shaw- Pynn  form  factor.  The  lengths 
of  the  vertical  lines  shown  here  are  proportional 
to  the  negative  of  the  weight  factors  for  an  intrin¬ 
sic  fault  in  an  (cc  crystal.  Figure  5  is  a  similar 
plot,  the  tabulated  I'„  in  this  case  being  those  of 
Animalu12  and  Yamamoto,11  while  the  correspond¬ 
ing  asymptotic  curve  was  computed  (or  the  Heine- 
Abarenkov  form  factor.  ,3~15  The  points  and  asymp¬ 
totic  curve  for  the  empty -core  model  correspond 
to  the  energy-wave- number  characteristic  shown 
in  Fig.  1.  Potential  functions  for  magnesium  and 
beryllium  were  also  computed  but  are  not  shown 
here  because  the  main  features  we  wish  to  con¬ 
sider  are  similar  for  all  three  polyvalent  metals. 

The  first  thing  we  should  note  regarding  Fig.  4 
is  that  it  is  very  easy  to  see,  from  a  plot  such  as 
this,  how  the  interactions  between  various  pairs 
of  planes  contribute  to  the  stacking-fault  energy. 

At  the  second  nearest  interplanar  distance  («  =  2) 
in  a'uminuin,  for  example,  J  is  positive  and  the 
weight  -N( 2)  is  also  positive.  This  indicates, 
according  to  Eqs.  (1)  and  (4),  that  the  interaction 
between  second- nearest -neighbor  planes  is  such 
as  to  oppose  the  formation  of  an  intrinsic  fault; 
i.  e. ,  this  particular  interaction  leads  to  a  positive 
contribution  to  the  stacking-fault  energy.  At  the 
next  interplanar  distance,  however,  the  contribu¬ 
tion  is  negative  because  J  and  -  A'(3)  are  of  oppo- 


FIG.  5.  Interplanar  potential  functions  for  aluminum 
based  on  the  Heim  -Abarenkov  and  empty -core-model  po¬ 
tentials. 


site  sign.  The  fourth-  and  fifth-neighbor  planes 
are  not  important  because  A'(4)  =  0  and  J(5h)*  0. 
Thus,  neglecting  possible  contributions  from  larg¬ 
er  n  (which  turns  out  to  be  valid  in  this  case),  we 
can  see  that  the  reason  aluminum  has  a  positive 
intrinsic  fault  energy,  or  equivalently,  the  reason 
that  the  fee  structure  is  stable  against  faulting,  is 
because  the  positive  contributions  of  second- near¬ 
est-neighbor  planes  outwe  gh  the  negative  con'rt- 
butions  of  third- nearest-  neighbor  planes.  The 
same  kind  of  analysis  can,  of  course,  be  applied 
to  the  magnesium  plot  shown  in  Fig.  3.  Here  the 
fact  that  the  hep  structure  is  stable  tells  us  tha* 
the  positive  stacking-fault-energy  contributions 
from  second-  and  thirc'-nearest-neighbor  planes 
must  be  greater  than  the  negative  contribution  front 
more  distant  pairs. 

Another  point  illustrated  in  Fig.  4  is  that  at  the 
interplanar  distances  of  importance  here,  it  makes 
very  little  difference  whether  we  use  the  data  of 
Shaw,  Shaw  and  Pynn,  or  Shyu  et  al.  in  the  calcu¬ 
lation  of  J .  Since  the  only  difference  among  these 
three  sets  of  data  is  the  way  in  which  exchange  and 
correlation  effects  are  handled  in  the  screening 
calculation,  the  comparison  showm  here  indicates 
that  such  effects  are  of  little  consequence  in  the 
present  context.  In  fact,  the  good  agreement  of  the 
results  obtained  with  the  other  two  characteristics 
with  the  result  based  on  Shaw's  data,  which  was 
calculated  in  the  Hartree  approximation,  suggests 
that  exchange  and  correlation  corrections  can  safe¬ 
ly  be  ignored  in  calculations  of  stacking- fault  en¬ 
ergies.  The  reason  is,  as  we  indicated  earlier, 
that  the  effects  of  such  corrections  are  observed 
at  interplanar  or  interatomic  distances  somewhat 
shorter  than  those  Involved  in  the  calculation  of 
stacking-fault  energies. 

Figure  4  also  illustrates  the  generally  good 
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agreement  experienced  between  the  numerically 
computed  value  of  J  and  the  value  obtained  Irom 
the  asymptotic  approximation  of  Blandine/  al.  at 
separation  distances  of  about  four  interplanar 
spacings.  Further  examples  are  provided  by  the 
Yamamoto  and  empty-core  results  shown  in  Fig. 

5.  Thus,  as  we  indicated  in  Sec.  Ill,  it  was  possi¬ 
ble  to  terminate  the  numerical- integration  com¬ 
putations  at  «  =  4  and,  in  effect,  use  the  asymptotic 
curve  for  larger  n  without  introducing  serious 
errors  in  most  applications. 

An  exception  is  illustrated  in  Fig.  5,  where  it 
is  evident  that  computed  values  based  on  Animalu’s 
semilocal  energy- wave- number  characteristic  are 
significantly  greater  than  the  corresponding  asymp¬ 
totic  curve.  This  difficulty  was  also  observed 
with  the  Animalu  calculation  for  magnesium  and, 
to  a  much  greater  extent,  with  the  beryllium  cal¬ 
culation.  In  fact,  with  beryllium,  the  disagree¬ 
ment  was  so  serious  that  we  were  not  able  to  ob¬ 
tain  even  a  reasonable  estimate  of  the  stacking- 
fault  energy  based  on  the  Animalu  data.  Some 
difficulty  was  also  experienced  with  the  beryllium 
potentials  based  on  the  data  of  Shaw,  Shaw-Pynn, 
and  Shyu  fl  al. ,  but  not  to  such  an  extent  as  to 
completely  invalidate  the  corresponding  estimates 
of  y. 

We  can  get  some  idea  of  the  influence  of  nonlo¬ 
cality  through  a  comparison  of  the  Animalu  and 
Yamamoto  results  shown  in  Fig.  5,  since  Yama¬ 
moto’s  data  were  obtained  from  a  fully  nonlocal 
treatment  of  the  same  model  potential  used  in 
Animalu's  calculations.  The  principal  result  is 
just  what  one  would  expect  from  Shaw's*  investi¬ 
gation  of  similar  effects  on  interionic  potentials, 
namely,  that  a  fully  nonlocal  treatment  tends  to 
diminish  the  amplitude  while  preserving  the  fre¬ 
quency  and  phase  of  the  long-range  oscillations  of 
the  potential  function.  It  should  be  noted,  however, 
that  Animalu’s  energy-wave-number  characteris¬ 
tic  is  semilocal,  in  the  sense  that  he  made  use  of 
certain  simplifying  approximations,  usually  invoked 
in  the  local  theory,  in  performing  the  final  integra¬ 
tion  to  obtain  F s.  12  If  a  truly  local  approximation 
had  been  employed,  Shaw’s*  studies  indicate  that 
the  amplitude  of  the  oscillations  in<7  would  be  even 
greater  than  that  obtained  with  the  Animalu  approxi¬ 
mation.  We  conclude,  therefore,  in  agreement 
with  Shaw,  that  nonlocality  must  be  accounted  for 
in  any  accurate  first- principles  predictions  of 
interionic  or  interplanar  potentials. 

This  is  not  to  say,  however,  that  we  cannot  use 
simpler,  empirically  adjusted  local  models  for  the 
prediction  of  trends  or  even  for  obtaining  rough 
estimates  of  the  values  of  the  stacking- fault  ener¬ 
gy  itself,  provided  that  we  have  some  prior  knowl¬ 
edge  of  the  general  features  of  the  nonlocal  energy- 
wave-number  characteristic.  This  is  illustrated 


by  the  fact  that  the  empty-core  calculations  shown 
in  Fig.  5  actually  agree  rather  well  with  the  more 
sophisticated  first-principles  results  shown  in 
Fig.  4.  The  energy-wave- number  characteristic 
used  here  was  of  the  usual  form,  11 

=  ll/€(?)jAfV'M.  (8) 

where  €(q)  is  the  Hartree  dielectric  function  and 
M(x)-  cos(irx/2x0)e'DxZ  .  (9) 

The  parameter  x0  was  fixed  by  requiring  that  the 
first  zero  in  F„  coincide  with  that  in  Animalu’s 
tabulation,  and  the  damping  parameter  D  was  cho¬ 
sen  to  bring  the  peak  near  q  -  2ft T  down  into  rough 
agreement  with  tabulated  nonlocal  functions.  The 
value  used  for  all  three  polyvalent  metals  was 
D  -  0. 12,  which  is  somewhat  larger  than  the  damp¬ 
ing  parameters  usually  employed  in  local-empiri¬ 
cal-pseudopotential  studies.  11  This  value  did,  how¬ 
ever,  lead  to  rather  good  agreement  with  the  non¬ 
local  energy  wavenumber  characteristics  and  in¬ 
terplanar  potential  functions  for  all  three  polyval¬ 
ent  metals,  the  results  shown  In  Figs.  1,  4  and 
5  being  typical. 

Turning  now  to  the  noble  metals  we  show  a  repre¬ 
sentative  result  in  Fig.  6.  This  is  the  interplanar 
potential  function  for  gold,  as  determined  from 
Moriarty's  energy- wave- number  characteristic,  IS 
along  with  the  corresponding  asymptotic  approxi¬ 
mation.  The  second  set  of  points  shown  here  was 
obtained  with  the  altered  energy-wave- number 
characteristic  shown  in  Fig.  2;  the  significance  of 
these  data  will  be  discussed  shortly.  Again,  the 
vertical  lines  are  proportional  to  the  weights 
-  N(n)  for  an  intrinsic  fault  in  an  fee  crystal. 

From  these  results,  and  those  for  copper  and 
silver  as  well,  it  is  evident  that  there  is  serious 
diagreement  between  the  interplanar  potential  func- 


FIG.  6.  Interplanar  poten'ia  functions  for  gold.  The 
corresponding  energy -wave-nu  obe-  characteristics  are 
shown  in  Fig.  2. 
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tion  predicted  by  the  asymptotic  theory  and  that 
obtained  by  direct  numerical  integration.  This  is 
not  surprising  in  view  of  the  fact  that  the  asymp¬ 
totic  approximation  is  based  on  the  assumption  that, 
as  far  as  the  long-range  behavior  of  J  is  concerned, 
the  logarithmic  singularity  at  q  =  2k r  is  the  domi¬ 
nant  feature  of  the  energy- wave- number  charac¬ 
teristic.  Thus,  when  the  asymptotic  form  resulting 
from  this  assumption  turns  out  to  be  a  rapidly  de¬ 
caying  exponential,  as  it  does  here,  there  emerges 
the  possibility  that  other  features  of  FM  may  domi¬ 
nate  at  large  values  of  nh.  Our  results  show  that 
this  is  in  fact  the  case  for  the  noble  metals  and 
suggest  that  similar  results  might  be  expected  for 
other  monovalent  metals  as  well. 

One  encouraging  consequence  of  the  radical  de¬ 
parture  of  the  numerically  computed  J  from  the 
asymptotic  approximation  is  that  the  computed 
values  lead  to  positive  energies  for  fee  faults  while 
the  asymptotic  theory  incorrectly  predicts  nega¬ 
tive  fault  energies. 5  As  is  evident  from  Fig.  6, 
this  happens  because  the  positive  contribution 
from  third- nearest- neighbor  planes  is  strong 
enough  to  offset  the  negative  contribution  from 
>i  =  2  interactions.  Thus,  the  strong  negative  dip 
in  J  between  «  =  2  and  3  seems  to  be  essential  to 
the  correct  prediction  that  the  fee  noble  metals  are 
stable  against  faulting. 

Unfortunately,  the  dip  in  J  is  not  strong  enough. 
In  Sec.  V  we  will  present  predictions  ol  y  that 
show  that  although  positive  energies  are  obtained 
for  fee  faults,  these  energies  are  low  by  about  an 
order  of  magnitude  for  gold  and  silver.  Further¬ 
more,  our  calculations  also  indicate  that  noble- 
metal  fault  energies  for  the  unstable  hep  crystals, 
which  one  would  expect  to  be  negative,  turn  out  to 
be  positive  and  quite  large.  This  indicates,  of 
course,  that  the  hep  structure  is  more  stable  than 
the  fee  structure,  an  incorrect  prediction  that 
agrees  with  Moriarty’s  total-energy  calculation. 1# 
The  trouble  is  obviously  due  to  some  inaccuracy 
in  the  Moriarty  F *  data,  probably,  as  he  suggests, 
to  his  neglect  of  the  subtle  effects  of  crystal  field 
splitting  on  s-d  hybridization.  16  However,  what¬ 
ever  the  basic  reason  may  be,  it  is  of  interest  to 
note  that  the  correct  interplanar  potential  function 
must  have  a  more  pronounced  minimum  than  that 
obtained  from  the  Moriarty  data.  This  would  pro¬ 
duce  a  greater  difference  between  the  positive 
(n  =  3)  and  the  negative  (n  =  2)  contributions  to  the 
fee  intrinsic-fault  energy,  while  at  the  same  time 
reducing  the  predicted  values  of  fault  energies  in 
the  hypothetical  hep  structure. 

An  Interplanar  potential  function  having  the  de¬ 
sired  characteristic  is  shown  in  Fig.  6.  This 
curve  was  generated  from  the  altered  F„  data 
shown  in  Fig.  2.  We  wish  to  emphasize  that  this 
alteration  is  entirely  arbitrary;  it  was  obtained  by 


trtal-and-error  manipulations  of  the  data  in  the 
neighborhood  of  q/kr=  x,  since  J  seems  to  be  most 
sensitive  to  the  detailed  curvature  of  F„  in  this 
region.  We  think  the  result  is  of  interest,  how¬ 
ever,  because  it  shows  that  for  the  noble  metals 
at  least,  seemingly  insignificant  changes  in  the 
energy-wave-number  characteristic  can  lead  to 
major  changes  in  the  interplanar  potential  (Figs.  2 
and  6)  and  in  the  predicted  stacking-fault  energy 
(the  alteration  shown  here  increased  the  fee  in¬ 
trinsic-fault  energy  by  about  a  factor  of  7).  Un¬ 
fortunately  this  extreme  sensitivity  to  detailed 
structure  in  the  energy- wave- number  characteris¬ 
tic  also  suggests  that,  in  contrast  to  the  situation 
for  polyvalent  metals,  there  is  little  hope  for 
making  useful  estimates  on  the  basis  of  simple 
models  for  the  noble  metals. 

V.  STACKING-FAULT  ENERGIES 
Table  n  is  a  listing  of  intrinsic-stacking-fault 
energies  for  the  stable  structures  of  aluminum, 
magnesium,  and  beryllium.  Calculations  were 
also  performed  for  extrinsic  and  twin  faults  and, 
with  few  exceptions,  the  approximate  relations  for 
hep  crystals10  yMl=  lyu«“  3y,.lB ,  were  found  to 
hold  within  about  10%.  The  exceptions  were  the 
twin-fault  energies  in  magnesium,  as  determined 
from  the  energy- wave- number  characteristics  of 
Shaw  and  Pynn,  Shaw,  and  Shyu  et  al, ,  where  we 


TABLE  11.  Intrinsic-stacking-fault  energies  in  erg/ 
cm2  for  aluminum  (fee),  magnesium  (hep),  and  beryllium 
(hep).  _  _ 


Data  source 

Al 

Mg 

Be 

Shaw— Pynn 

57.  9 

13.  9 

439 

Shaw 

52.  1 

8.  54 

390 

Shyu  et  al. 

69.2 

16.6 

468 

Animalu 

104 

33.4 

... 

Yamamoto 

110 

37.7 

305 

Empty-core  model 

41).  1 

11.5 

128 

Other  calculations 

195* 

62b 

142' 

8.7f 

30* 

50* 

225* 

Experimental 

estimates 

135d 

135—280* 

60— 93h 
100-1 50d 

190' 

•Reference  17;  based  on  Animalr  data. 

•Reference  18;  based  on  Shaw  data. 

‘Reference  18;  based  on  Aw  mala  data. 

•Reference  20. 

•P.  C.  J.  Gallagher,  Met.  Trans.  1,  2429  (1970). 
'Reference  19;  based  on  Shav.  data. 

•Reference  9;  based  on  Hamsun's  pseudopotential. 
*D,  H.  Sastry,  Y.  V.  R.  K.  ’rasad,  and  K.  1.  Vasu. 
Scripts  Met.  9,  927  (1969). 

'H.  Conrad,  Air  Force  Mater  air  Laboratory  Report 
No.  AFML-1R-65-310  (1965)  (unpublished). 
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TABLE  III.  Intrinsic-stAcking-fault  energies  in  erg/ 
cm1  (or  fee  copper,  silver,  and  gold.  Calculated  values 
are  based  on  the  energy-wave-number  characteristics 
of  Moriarty  (Ref.  16). 


Cu 

Ag 

Au 

Calculated 

44.8 

1.7 

4.8 

Average 

experimental 

estimate* 

55 

21.7 

50 

•P.  C.  J.  Gallagher,  Met.  Trans.  1,  2429  (1970). 


found  y„t  =  2.  3ytwl>.  With  fee  aluminum  we  ob¬ 
tained,  to  within  about  2Sfc,  y.,t  =  yi.t  =  2y,.u .  For 
the  unstable  close- packed  structures  (hep  alumi¬ 
num  and  fee  magnesium  and  beryllium)  we  obtained 
negative  energies  for  all  three  types  of  fault,  thus 
correctly  indicating  the  instability  of  these  struc¬ 
tures  against  faulting. 

Intrinsic-fault  energies  for  fee  copper,  silver, 
and  gold,  as  determined  from  Moriarty’s  energy- 
wave-number  data,  are  shown  in  Table  IIL  As 
was  noted  in  the  Sec.  IV,  the  hep  fault  energies 
for  these  metals,  rather  than  being  negative  as 
they  should  be,  were  calculated  to  be  positive  and 
quite  large  (y-  26,  15,  and  25  erg/cm2  for  the  hep 
intrinsic-fault  energies  in  copper,  silver,  and 
gold,  respectively).  The  good  agreement  with  ex¬ 
periment  obtained  for  the  fee  intrinsic-fault  ener¬ 
gy  in  copper  is  therefore  probably  fortuitous. 

When  the  energy- wave- number  characteristic 
for  gold  was  altered  as  shown  in  Fig.  2,  we  ob¬ 
tained  37  and  20  erg/cm2  for  the  fee  and  hep  intrin¬ 
sic  faults,  respectively.  We  wish  to  emphasize 
again,  however,  that  although  this  change  is  in  the 
right  direction,  the  alteration  of  F„  shown  in  Fig. 

2  is  entirely  arbitrary,  and  serves  only  to  illus¬ 
trate  the  sensitivity  of  noble- metal  results  to 
rather  minor  features  in  the  energy- wave- number 
characteristics.  This  being  the  case,  there  is 
not  much  more  we  can  say  about  the  noble  metals 
except  to  note  that  the  calculation  of  stacking-fault 
energies  appears  to  be  a  very  demanding  test  of 
the  accuracy  of  the  energy- wave- number  charac¬ 
teristic  in  the  vicinity  of  the  smallest  reciprocal- 
lattice  spacing. 

Returning  now  to  the  polyvalent  metals  we  can 
see,  as  anticipated  from  the  plots  of  the  inter- 
planar  potential  function  shown  in  Fig.  4,  that  ex¬ 
change  and  correlation  effects  are  not  particularly 
significant  since  little  change  was  experienced  in 
going  from  the  Hartree  approximation  (Shaw)  to 
calculations  with  such  many-body  corrections  in¬ 
cluded  (Shaw  and  Pynn  and  Shyu  el  al. ).  The  only 
appreciable  effect  was  that  shown  for  magnesium, 
where  there  is  roughly  a  factor  of  2  spread  in  pre¬ 
dicted  intrinsic-fault  energies. 


Results  based  on  the  energy- wave- number  data 
of  Animalu  and  Yamamoto  indicate  that,  in  two 
cases  at  least,  differences  between  fully  nonlocal 
and  semilocal  calculations  are  not  as  great  as  one 
might  expect  from  the  data  shown  in  Fig.  6.  Ap¬ 
parently  what  has  happened  in  aluminum  is  that 
the  large  differences  between  the  Animalu  and 
Yamamoto  calculation  of  J  at  second-  and  third- 
nearest-  neighbor  planes  are  approximately  can¬ 
celed  in  performing  the  sum  in  Eq.  (7)  (recall  that 
the  second- nearest -neighbor  plane  gives  a  positive 
contribution  to  y  while  the  third-nearest- neighbor 
plane  yields  a  negative  contribution).  A  similar 
situation  exists  for  magnesium,  the  cancellation 
in  this  case  being  between  the  positive  contribu¬ 
tions  of  the  second-  and  third-neighbor  pairs  and 
the  negative  contribution  of  fourth- nearest- neigh¬ 
bor  planes.  For  beryllium,  however,  the  effect 
was  quite  marked  because  the  Animalu  data  gave 
correction  terms  [those  in  the  sum  on  the  right 
side  of  Eq.  (7)]  that  were  about  an  order  of  magni¬ 
tude  greater  than  the  corresponding  values  obtained 
with  Yamamoto’s  data.  As  was  mentioned  pre¬ 
viously,  it  is  this  large  deviation  of  the  numeri¬ 
cally  integrated  J  from  the  asymptotic  approxima¬ 
tion  that  prevented  us  from  obtaining  a  meaningful 
estimate  of  y  for  beryllium  from  Ammalu’s  data. 
One  should,  of  course,  expect  to  encounter  diffi¬ 
culties  with  the  local  or  semilocal  approximation 
for  beryllium  because  the  absence  of  />  states  in 
the  core  leads  to  a  strongly  nonlocal  pseudopoten¬ 
tial.  11 

It  would  appear,  however,  from  results  such  as 
those  illustrated  in  Figs.  1,  4,  and  5,  that  in  the 
present  context  the  principal  effect  of  a  fully  non¬ 
local  treatment  is  the  suppression  of  the  peak  in 
F „  near  q-  Zkr  and  the  resulting  reduction  in  the 
interplanar  potential  function  <7.  Thus,  given  some 
prior  knowledge  of  the  magnitude  of  the  peak  in 
the  nonlocal  F„ ,  it  should  be  possible  to  use  em¬ 
pirically  adjusted  local- model  potentials,  even  for 
such  strongly  nonlocal  elements  as  beryllium,  to 
obtain  rough  estimates  of  y.  To  test  this  idea,  we 
used  the  modified  empty-core  model  described  in 
the  Sec.  IV  to  compute  the  intrinsic-fault  energies 
shown  in  Table  II.  Considering  the  fact  that  the 
model  is  extremely  simple  and  that  the  same  damp¬ 
ing-factor  adjustment  was  used  for  all  three  met¬ 
als,  the  agreement  with  other  predictions  of  v  is 
surprisingly  good.  Although  improved  agree¬ 
ment  could  probably  be  obtained  by  further  rdjust- 
ments  of  the  model-potential  parameters,  we  see 
little  reason  for  doing  so,  particularly  since  there 
is  still  considerable  uncertainty  as  to  what  one 
should  assume  for  the  correct  value  of  v.  In  any 
case,  even  if  y  were  known  more  accurately,  the 
real  merit  of  such  simplified  models  lies  in  their 
ability  to  predict  trends  in  relative  magnitudes  •'•■d. 
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perhaps,  obtain  rough  estimates  ot  y  itself,  with 
minimum  requirement  for  prior  adjustment  of  the 
model-potential  parameters.  In  this  sense,  then, 
the  rough  agreement  shown  here  is  actually  quite 
encouraging. 

As  far  as  the  more  rigorous  models  are  con¬ 
cerned,  comparisons  of  our  computations  with 
those  performed  by  other  methods  met  with  mixed 
results.  With  aluminum,  for  example,  our  value 
of  104  erg/cm2  based  on  Animalu’s  data  is  signifi¬ 
cantly  lower  than  Hodges17  value  of  195  erg/cm1, 
which  was  computed  from  the  same  energy-wave- 
number  characteristic  using  the  reciprocal- space 
formulation.  It  is  also  lower  than  the  value  of 
142  erg/cm2  obtained  by  Wilkes  and  Sargent,  “  who 
used  the  same  data  in  a  real- space  treatment  in¬ 
volving  the  sums  of  inter  ionic  potentials  over  a 
large  number  of  pairs  of  ions  in  the  faulted  and 
perfect  structures.  There  is,  unfortunately,  rea¬ 
son  to  question  the  accuracy  of  both  the  Wilkes- 
Sargent  result  and  our  own.  The  Wilkes-Sargent 
calculation  did  not  make  use  of  the  asymptotic 
theory  of  Blandinet  al.  and,  therefore,  suffers 
from  a  convergence  difficulty  typical  of  pairwise 
sums  over  the  direct  lattice.  Our  calculation  in 
this  case  is  also  subject  to  a  convergence  uncer¬ 
tainty  because,  as  we  noted  earlier,  with  the 
Animalu  data  it  was  not  possible  to  carry  the  cal¬ 
culation  out  to  distances  large  enough  to  obtain 
agreement  with  the  asymptotic  potential.  We  can 
be  less  certain  of  possible  sources  of  error  in 
Hodges  result,  although  the  fact  that  he  had  to 
perform  a  number  of  rather  sensitive  computations 
by  hand  could  have  resulted  in  some  inaccuracy. 
The  only  other  comparison  we  can  make  for  alu¬ 
minum  is  our  value  of  52  erg/cm2,  obtained  from 
Shaw’s  data,  with  the  Wilkes-Sargent  value  of  62 
erg/cm2.  In  view  of  the  convergence  difficulty 
experienced  by  Wilkes  and  Sargent  this  ageement 
is  probably  as  good  as  can  be  expected. 

For  magnesium,  our  results  based  on  the  Shaw 
and  Animalu  energy-wave-number  data  are  in 
excellent  agreement  with  the  values  reported  by 
Ducharme1’  (8.  7  erg/cm2  with  Shaw’s  data)  and 
Hodges17  (30  erg/cm2  with  Animalu’s  data).  Har¬ 
rison’s  calculation  of  50  erg/cm2  was  based  on  a 
different  energy- wave- number  characteristic’  and 
is  included  here  only  as  a  further  illustration  of  the 
rather  significant  differences  that  can  result  from 
the  use  of  different  potentials.  The  only  compari¬ 
son  we  can  make  for  beryllium  is  our  calculation 
of  the  fault  energy  with  that  of  Ducharme, 12  both  of 
which  made  use  of  Shaw’s  energy- wave- number 
characteristics.  Here  again,  however,  we  ex¬ 
perienced  some  convergence  difficulty,  which  could 
account  for  the  discrepancy. 

Agreement  with  experimental  results  is  general¬ 
ly  rather  poor.  We  are  not  particularly  concerned 


that  our  results  for  beryllium  seem  to  be  too  high, 
because  of  the  convergence  problem  discussed 
earlier.  The  poor  agreement  for  magnesium  and 
aluminum,  however,  deserves  further  comment. 

First  we  should  point  out  that  it  is  quite  possible 
that  the  experimental  values  for  these  metals  are 
too  high,  since  experimental  difficulties  are  usually 
encountered  when  y  is  greater  than  about  20  erg/ 
cm2.  20  The  source  of  this  difficulty  is  that  fault 
widths  are  roughly  proportional  to  the  inverse  of 
the  stacking-fault  energy  and  are  therefore  too 
small  for  direct  observation  when  y  is  large.  One 
must  then  resort  to  indirect  methods,  such  as  the 
observation  of  dislocation- loop  annealing  rates, 20 
to  obtain  estimates  of  y.  Interpretations  of  such 
data  necessarily  involve  the  assumption  of  physical 
models  and  subsequent  calculations  to  extract  y 
from  the  direct  experimental  results.  The  pres¬ 
ent  situation  is  therefore  one  of  considerable  un¬ 
certainty,  as  evidenced  by  the  spread  in  experi¬ 
mental  data  referenced  in  Table  I.  For  this  rea¬ 
son  we  are  inclined  to  regard  comparisons  with 
experiment  primarily  as  tests  of  the  ability  to  pre¬ 
dict  trends,  placing  less  emphasis  on  the  achieve¬ 
ment  of  absolute  numerical  agreement. 

As  for  why  the  calculations  might  be  inadequate, 
if  indeed  it  is  the  calculations  that  are  at  fault, 
we  cannot  offer  a  completely  satisfactory  explana¬ 
tion.  While  it  is  certainly  possible  that  the  nu¬ 
merical  problems  discussed  in  Sec.  II  led  to  sig¬ 
nificant  errors  in  the  predictions,  agreement  with 
calculations  performed  by  other  methods  seems 
to  indicate  that  this  is  not  the  case.  Nor  does  it 
seem  likely  that  uncertainties  in  exchange  ard 
correlation  corrections  are  responsible,  because 
we  can  totally  ignore  such  effects  and  still  ob’air 
about  the  same  answers.  Nonlocal  effects  can  be 
ruled  out  because  they  are  fully  accounted  for  in 
some  of  the  energy- wave- number  characteristics 
of  concern  here.  Although  magnesium  and  alumi¬ 
num  are  generally  considered  "simple”  metals, 
in  the  sense  that  the  matrix  elements  of  their 
pseudopotentials  are  small,  we  cannot  ignore  the 
possibility  that  higher-order  terms  in  the  pertur¬ 
bation  expansion,  which  would  be  structure  depen¬ 
dent  and  therefore  different  in  the  faulted  and  per¬ 
fect  crystals,  may  have  a  greater  effect  in  stack¬ 
ing-fault-energy  calculations  than  they  do  in  cal¬ 
culations  of  other  properties.  Still,  if  agreement 
with  presently  available  experimental  data  is  our 
goal,  we  are  looking  for  a  factor  of  about  2  or  3 
increase  in  y,  and  the  possibility  that  a  higher- 
order  perturbation  treatment  would  give  us  such  a 
factor  does  seem  remote. 

Thus,  because  the  major  sources  of  error  usually 
present  in  defect  studies  are  of  diminished  impor¬ 
tance  in  stacking-fault-energy  calculations,  it  is 
difficult  to  identify  any  one  uncertainty  as  being 
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most  likely  to  cause  disagreement  with  experi¬ 
ment.  The  best  we  can  do  at  present  is  to  simply 
note  that  there  are  still  many  approximations  and 
idealizations  present,  both  in  the  theory  and  in  the 
derivation  of  stacking-fault  energies  from  experi¬ 
mental  observations  and,  for  this  reason,  perhaps 
we  should  not  expect  much  better  agreement  than 
that  obtained  here. 

VI.  SUMMARY  AND  CONCLUSIONS 

We  have  applied  the  method  of  Blandin,  Friedel, 
and  Saada*  to  the  calculation  of  stacking-fault 
energies  in  aluminum,  magnesium,  beryllium, 
copper,  silver,  and  gold.  Energy-wave-number 
characteristics  representing  various  degrees 
of  refinement  of  the  theory  were  used  in  order 
that  the  sensitivity  of  the  result  to  such  refine¬ 
ments  could  be  determined. 

Our  results  for  the  polyvalent  metals  indicate 
that  the  stacking-fault  energy  is  influenced  by 
many-body  and  nonlocal  effects,  but  not  nearly  so 
much  as  to  invalidate  the  use  of  approximate 
model  potentials.  In  fact,  comparisons  of  results 
based  on  a  simple  local  potential  with  those  based 
on  the  more  sophisticated  nonlocal  treatments  sug¬ 
gest  that,  with  certain  semiempirical  adjustments 
to  account  for  the  major  effects  of  nonlocality,  the 
simpler  theory  is  entirely  adequate  for  predicting 
trends  and  for  obtaining  rough  estimates  of  the 
stacking-fault  energy,  as  long  as  the  use  of  such 
models  is  restricted  to  polyvalent  metals.  Cal¬ 
culations  for  copper,  silver,  and  gold,  on  the 
other  hand,  show  that  the  stacking-fault  energies 
of  these  metals  (and  probably  other  monovalent 
metals  as  well)  are  quite  sensitive  to  minor  altera¬ 
tions  in  the  energy-wave-number  characteristic. 
The  use  of  simplified  semiempirical  models  for 
these  metals  could  therefore  lead  to  serious  error. 

The  real- space  formulation  of  Blandin  et  al. 
was  found  to  have  both  its  advantages  and  dis¬ 
advantages.  On  the  positive  side,  the  fact  that  the 
method  involves  a  sum  of  interplanar  potentials 
over  pairs  of  close-packed  planes  offers  some 
interpretative  advantages  over  the  reciprocal- 
space  treatment.  It  is  very  easy  to  see,  for  ex¬ 
ample,  in  terms  of  the  amplitudes  and  oscillatory 
characteristics  of  the  interplanar  potentials,  how 
interactions  between  various  pairs  of  planes  con¬ 
tribute  *o  the  stability  or  instability  against  faulting 


for  a  given  close- packed  structure. 

The  principal  disadvantage  of  the  method  is  that 
it  presents  some  difficult  computational  problems, 
both  in  the  sensitivity  of  the  result  to  numerical- 
integration  approximations  and  in  the  usual  con¬ 
vergence  problem  associated  with  real-space  lat¬ 
tice  sums.  Still,  with  the  aid  of  the  asymptotic 
approximations  developed  by  Blandin  et  al. ,  the 
latter  difficulty  can  be  largely  overcome,  and  the 
method  therefore  offers  an  instructive,  though 
perhaps  somewhat  less  accurate.  rnative  to 
the  more  familiar  reciprocal-sp..  atmont. 

Comparisons  of  predicted  stac!  -  ’  rgies 

with  experimental  results  were  disa  iiw 
calculated  values  in  most  cases  bet  :  ^ubstant: 
ly  lower  than  experimental  values,  o.nc.  c 
major  uncertainties  that  one  usually  has  v 
with  in  defect-energy  calculations  (the  v  . 

perturbation  theory,  exchange  and  ccr: i at ■ 
corrections,  and  the  importance  of  no nloca-  >  ts) 
are  of  minimal  significance  in  the  present  appli¬ 
cation,  it  is  difficult  to  explain  why  this  disagree¬ 
ment  exists.  Although  this  is  an  unsettling  situa¬ 
tion,  and  one  that  we  think  deserves  further  study, 
we  can  draw  some  satisfaction  from  the  fact  that 
the  calculations  do  correctly  predict  the  stabilities 
of  close- packed  structures  against  faulting,  and, 
at  least  roughly,  the  relative  magnitudes  of  stack¬ 
ing-fault  energies  for  different  metals.  The  cal¬ 
culations  also  indicate  that  trends  such  as  these 
can  be  predicted  with  reasonable  success  with  a 
local- model  potential,  provided  that  the  major 
effects  of  nonlocality,  the  suppression  of  the  peak 
in  the  energy- wave- number  characteristics  at  q 
-  2kr,  is  approximately  accounted  for  by  an  em¬ 
pirical  adjustment.  This  may  well  be  the  most 
important  result  of  the  present  work  because  it 
suggests  that  in  spite  of  quantitative  difficulties 
that  seem  to  persist  in  even  the  most  refined  ver¬ 
sion  of  the  theory,  there  is  still  much  that  can  be 
learned  from  the  simplest  local  approximation. 
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An  energy-dependent  f -matrix  generalization  of  Harrison's  theory  of  multi-ion  interactions  is 
developed.  Explicit  expressions  are  derived  for  two-  and  three-ion  potentials  in  a  system  of  ions  having 
{/-electron  scattering  resonances  of  the  Breit-Wigner  form.  It  is  shown  that  although  a  simple  pairwise 
interionic  potential  does  not  exist,  one  can  define  an  effective  potential  such  that  the  correct 
band-structure  energy  is  obtained  in  the  form  of  a  pairwise  sum.  The  effective  pairwise  potential  is  then 
computed  through  third  order  for  an  idealized  transition-metal  senes.  The  results  show  that  third-  and 
possibly  higher-order  interactions  must  be  included  in  the  calculation  to  fully  account  for  the  effects  of 
d  -band  splitting  on  interionic  forces. 


I  INTRODUCTION 

The  underlying  assumption  in  the  pseudopotential 
theory  of  nontransition  metals  is  that  conduction  - 
electron-ion  scattering  is  weak  enough  to  be  treated 
in  the  first  Born  approximation.1,2  On  the  basis  of 
this  assumption,  one  then  calculates  electro,  den¬ 
sities  and  screening  potentials  that  are  self-consis¬ 
tent  to  first  order  in  the  electron -ion  pseudopoten¬ 
tial  leading,  via  standard  perturbation  theory,  to 
an  expression  for  the  total  electronic  energy  that 
is  correct  to  second  order  in  the  pseudopotential. 

By  taking  the  Fourier  transform  of  the  second-or¬ 
der  term  in  this  expression,  one  then  obtains  a 
rather  simple  formula  for  the  electronic  part  of  the 
interionic  potential.  As  it  turns  out,  at  this  level 
of  approximation  the  resulting  potential  is  indepen¬ 
dent  of  the  positions  of  neighboring  ions;  the  only 
geometrical  factor  that  enters  is  the  separation 
distance  of  the  ion  pair. 

The  application  of  this  simple  model  to  transition 
metals  is,  however,  clearly  invalid.  Here  it  is 
well  known  that  the  indirect  interionic  interaction 
is  dominated  by  resonant  scattering  of  d  electrons 
and  that  a  low -order  perturbation  treatment  is 
therefore  quite  unrealistic.  One  might  hope  to  im¬ 
prove  the  calculation,  while  still  retaining  the  for¬ 
mal  simplicity  of  the  standard  theory,  by  starting 
with  an  exact  calculation  for  scattering  by  individu¬ 
al  ions,  while  keeping  the  second-order  approxi¬ 
mation  to  the  contributions  of  interionic  scattering 
to  the  total  energy.  However,  as  Heine  and  Weaire 
have  noted, 2  this  leads  to  an  interaction  energy  that 
is  correct  only  for  an  isolated  pair  of  ions;  the 
convenient  pairwise  interaction  form  is  thus  ob¬ 
tained  by  simply  ignoring  the  rest  of  the  crystal. 

In  view  of  the  fact  that  scattering  by  neighboring 
ions  is  known  to  have  a  profound  effect  on  even  the 
atomiclike  aspects  of  the  electronic  structures  of 
transition  metals  (e.  g. ,  the  splitting  of  a  single¬ 
ion  d  resonance  into  subbands1),  it  would  seem  that 
the  isolated-pair  model  could  lead  to  serious  error 
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when  applied  to  the  calculation  of  interionic  inter¬ 
actions  in  real  crystals.  What  is  needed,  there¬ 
fore,  is  a  way  of  going  beyond  second  order  in  in¬ 
terionic  scattering  in  order  to  properly  account  for 
the  effects  of  neighboring  ions  on  interionic  inter¬ 
actions. 

Our  main  concern  in  the  present  work  is  the  de¬ 
velopment  of  a  formalism  appropriate  to  such  high¬ 
er-order  calculations  for  transition  metals.  The 
approach  we  use  is  similar  to  that  used  by  Harri¬ 
son  in  the  development  of  his  theory  of  multi -ion 
interactions  in  nontransition  metals.4  The  basic 
idea  is  that  one  can  write  the  perturbation  expan¬ 
sion  of  the  structure -dependent  part  of  the  elec¬ 
tronic  energy  in  a  form  such  that  the  leading  term 
is  the  sum  of  pairwise  interionic  potentials  and 
successive  terms  are  sums  of  three -ion  and  high¬ 
er-order  interactions.  Each  multi -ion  interaction 
term  in  this  expansion  corresponds  to  a  closed 
scattering  path  connecting  a  group  of  ions  and  the 
expansion  is  therefore  conveniently  represented  as 
a  sum  of  multi-ion  interaction  diagrams.  Within 
this  framework  the  higher -order  correction  terms 
we  seek  are  therefore  those  corresponding  to  dia¬ 
grams  with  three  or  more  vertices. 

There  are,  however,  two  generalizations  of  Har¬ 
rison’s  formalism  that  are  required  before  the 
multi -ion  expansion  idea  can  be  applied  to  transi¬ 
tion  metals.  The  original  development,  because 
it  was  intended  for  applications  to  nontransition 
metals,  was  based  on  the  assumption  of  weak,  en¬ 
ergy  -independent  scattering  by  individual  ions. 

Both  approximations  are,  of  course,  grossly  inac¬ 
curate  in  transition  metals  where  near -resonance 
conditions  at  energies  of  interest  make  scattering 
quite  strong  and  highly  energy  dependent. 

The  extension  to  a  partial  wave  treatment  of 
scattering  by  individual  ions,  which  permits  one 
to  avoid  the  weak  scattering  assumption,  is  actual¬ 
ly  quite  straightforward.  Such  a  development  was 
reported  first  by  Harrison  in  a  note  added  to  his 
original  article,4  and  later  by  Lloyd  and  Oglesby,* 
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who  based  their  development  on  a  different  expan¬ 
sion  technique.  An  alternative  treatment,  which 
more  closely  resembles  Harrison’s  original  ap¬ 
proach,  is  contained  in  Sec.  II  of  this  paper. 

It  is,  however,  the  treatment  of  energy -depen¬ 
dent  scattering  that  presents  some  difficulty  and 
is,  therefore,  of  more  concern  here.  Actually, 
there  are  two  parts  to  this  problem,  the  first  being 
encountered  in  the  early  stages  of  the  formal  de¬ 
velopment  and  the  second  in  handling  an  integral 
over  a  product  of  resonant  scattering  amplitudes. 

If  one  starts  with  an  expansion  similar  to  Harri¬ 
son’s,  the  formal  difficulty  occurs  when  the  energy- 
independent  pseudopotential  form  factor  used  by 
Harrison  is  replaced  by  energy -dependent  elements 
of  the  corresponding  /  matrix.  In  Sec.  II  we  show 
that  this  can  be  handled  by  combining  certain 
classes  of  scattering  diagrams.  With  the  Lloyd- 
Oglesby  expansion,  incidentally,  this  formal  dif¬ 
ficulty  is  avoided  altogether,  but  at  the  expense, 
we  think,  of  some  loss  of  clarity  as  to  which  dia¬ 
grams  are  to  be  included  in  the  sum.  Aside  from 
this  point,  however,  our  formal  result  is  in  exact 
agreement  with  that  of  Lloyd  and  Oglesby. 

In  Sec.  Ill  we  address  the  second  part  of  the 
problem,  that  of  performing  an  integral  over  ener¬ 
gy-dependent  scattering  amplitudes.  To  do  this  we 
separate  the  resonant  and  nonresonant  parts,  as¬ 
sume  the  Breit-Wigner  form  for  the  resonant  part, 
and  make  use  of  some  of  the  approximations  em¬ 
ployed  by  Mezei  and  Grime  r®  in  the  treatment  of  a 
similar  problem  to  perform  the  resulting  integrals 
analytically.  The  end  result  is  a  general  expres¬ 
sion  for  the  asymptotic  form  of  the  multi-ion  inter¬ 
action  energy  for  cases  where  ion  scattering  am¬ 
plitudes  contain  both  resonant  and  nonresonant 
parts. 

Next  we  turn  to  the  problem  posed  earlier  in  this 
discussion,  namely,  an  investigation  of  the  role 
played  by  multiple  interionic  scattering  in  deter¬ 
mining  interionic  interactions  in  a  crystalline  en¬ 
vironment.  First,  we  show  that  one  can  define  an 
effective  pairwise  potential  that  includes  multi -ion 
contributions.  This  effective  interionic  potential, 
which  necessarily  depends  on  crystal  structure  and 
on  the  orientation  of  the  ion  pair  with  respect  to 
crystallographic  axes,  is  shown  to  consist  of  the 
usual  potential  for  an  isolated  pair  plus  correction 
terms  corresponding  to  scattering  by  neighboring 
ions.  It  therefore  forms  a  suitable  basis  for  the 
study  of  interionic  forces  in  real  crystals. 

Finally,  in  Sec.  IV,  we  present  some  numerical 
computations  of  the  effective  interionic  potential 
for  nearest  neighbors  in  a  hypothetical  face-cen- 
tered-cubic  transition -metal  series.  The  results 
show  that  as  the  Fermi  level  is  varied  from  the 
bottom  to  the  top  of  the  d  band,  the  interaction  en¬ 
ergy  for  an  isolated  pair  goes  through  a  strong 


negative  peak  in  the  vicinity  of  the  resonance  ener¬ 
gy  followed  by  a  weaker  repulsive  peak  at  higher 
energies.  The  interesting  point,  however,  is  that 
when  third -order  interactions  with  neighboring  ions 
are  added,  the  structure  of  the  attractive  peak  is 
modified  in  such  a  way  as  to  indicate  the  presence 
of  an  additional  repulsive  component  at  energies 
just  above  resonance.  This  is  due,  we  suggest, 
to  changes  in  the  structures  of  bonding  and  anti- 
bonding  subbands  brought  about  by  multi-ion  scat¬ 
tering.  Thus,  as  we  indicated  above,  it  would 
seem  that  third-  and  higher-order  scattering 
among  neighboring  ions  does  indeed  play  a  signifi¬ 
cant  role  in  determining  the  nature  of  interionic 
interactions  in  transition -metal  crystals. 

II.  FORMAL  THEORY 

Our  starting  point  is  Harrison’s  expression  for 
the  configuration-dependent  part  of  the  band -struc¬ 
ture  energy4 

V=-  fBr  /'  p(E')dE'dE  ,  (1) 

Jo  Jo 

where  p(£)  is  the  density  of  states, 

p(£)  =  -  (2/ ir)Im  TrG(£j  . 

Here  Tr  denotes  the  trace  and  G(E)  is  the  Green’s 
operator  for  electrons.  From  the  well-known  re¬ 
lation 

G(E)  =  G0(E)  +  G0(E)T(E)G0(E)  ,  (2) 

we  see  that  t3  contains,  in  addition  to  the  configura¬ 
tion-dependent  term  containing  the  crystal  /  ma¬ 
trix  7(E),  a  free-electron  contribution  that  is  in¬ 
dependent  of  the  arrangement  of  ions  (G0  is  the 
free-electron  propagator).  Since  we  are  interested 
in  structure -dependent  properties,  we  will  drop 
this  term  and  calculate  only  the  contribution  of  the 
second  term  in  Eq.  (2). 

Next  we  use  the  /-matrix  expansion7 

r=E  T„ ,  (3) 

with 

r„=E  E  •••  E  (4) 

Rj  Hj<R1 

where  /(R,)  is  the  /  operator  for  an  ion  at  Rr  But 
the  first  term  in  Eq.  (3)  is  also  independent  of  ion 
arrangement  and  can  therefore  be  omitted.  This 
leaves,  for  the  structure -dependent  energy, 

V=  -  -E  [*'  f*  ImTrG07-nG0d£'rf£:  .  (5) 

*  n*  Jo  '  0 

By  making  use  of  the  invariance  of  the  trace  un¬ 
der  cyclic  permutation  of  the  operators,  and  from 
the  identity 


R  .  E  . 


BE  ISSNER 


.n  [with  //  =  /(S<)] 

'o7„Go  =  E  E  •••  E  TrG0tiC0t2  •  •  •  G0t„G0  =  -E  E  *••  E  Tr  — ®  (iG0tz  •  ■  •  G0  tn 

Bl  R1>>R2  Rn«n-1  Si  #2*51  Rn<R„.i 

R„*Ri 

-E  E  •••_  E  TrCo^!  tiG0 /2 * •  •  Go t^-i  •  ^ 

Bi  a2*8i  5n_i*B„^ 

R«-l*Rl 


;t  term  on  the  right-hand  side  is  the  sum 
n -ion  scattering  paths  with  the  first  and 
s  on  different  sites  while  the  second  term 
onds  to  diagrams  with  the  first  and  last  ions 
■ame  site.  The  next  step  is  to  show  that  this 
class  of  diagrams  can  be  combined  with 
irms  to  express  the  sum  in  Eq.  (5)  in  a 
onvenient  form . 

let  v  denote  the  energy-independent  poten- 
a  particular  ion,  then  the  corresponding  t 
ir  satisfies 


')=(  +  iG0(E')/(£') 


(7) 


( - - - — 

from  which  we  obtain 

‘%P=[l-vC0{S'))-'  *.§?/(£');  (8) 

but  from  Eq.  (7)  we  have 
[l-t-GolE')]-1*  =  t(E') 
and  Eq.  (8)  can  therefore  be  written 


dt 


(9) 


Next  we  substitute  Eq.  (9)  in  Eq.  (6)  and  perform 
cyclic  permutations  of  indices  and  operators  to 
show  that 


J 


G0T„G0= -E  E  •••  E  TrG0/iG0  •  •  •  /*  •  •  •  G0/„ 

B,  Bjrti  VS,.! 

Rn*Rl 

-  ^  "  "  a  ^  TrG0  t\C9  G0  •  •  •  C0  tx  , 

Ri  Ra'Ri  s„-i's„-2  d£ 

Bn-lrfR  t 

*  is  any  index  in  the  appropriate  range.  Since  there  are  n  such  equal  expressions  for  the  first  term 
- 1  expressions  for  the  second  term,  we  can  add  and  divide  by  the  numbers  of  expressions  to  obtain 

G0T„G0  =  TrG0(T^  +  Tna>)G0  ,  (11) 


(10) 


rG0T^Gt=-lnE  E^...Ejv(^  fJG0...^  +  G0<1g?t2...tB  +  ...+G,fl...gRfn)  (12) 


Ri»-1*R1 


olds  for  all  n  >  2.  The  pairwise  interaction 
n  -  2 )  is 


(13) 


Tr  -jjp-  (G0  itG0 12) 
*  B!  fijrf,  “ 


jGo 


This  process  can  be  continued  for  all  higher-order 
terms,  with  the  result 


tiGo  tg  +  C8  h 


dCy 

dE' 


m 

E 

r-2 


TrG0rsG0 


an  be  combined  with  Eq.  (13)  for  n  =  3  to  give 
V7?’)G0 


-E^EE  • 


(G©/  j  G8/2  • . .  Go/„). 

(14) 
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From  this  point  on  the  development  is,  for  the  most  part,  the  same  as  Harrison’s  treatment.  We  will, 
therefore,  give  only  an  abbreviated  sketch,  concentrating  on  those  points  where  differences  do  occur. 

To  calculate  the  trace  in  Eq.  (14)  we  express  all  operators  in  terms  of  their  plane -wave  matrix  ele¬ 
ments  .  Thus 


Trtoo^CVj •  •  •  GqO'E  E  •  •  •  E 

‘1  ‘2  ‘n 


i  t\  i  £g)  ( £2 1  tj  I  ) » *  ■  ( [  tn  I  Kj ) 
<£*  -  *?)(£’  -  £)  •••(£'- 


(15) 


where  b"  --  E'  +  ;5  and  (x  I  C>  =  O'1  /2  e1*"1  with  O  being 
the  volume  of  the  crystal.  Since  t s  is  the  t  operator 
for  an  ion  at  R,  we  have 

<£,  !  I ,  S,.l)  =  eni^,-R/<£y  j 

=  _  e“5y.>-S 2m)  , 

(16) 

— - — — _ l 


r - — 

where  tg  is  the  t  operator  for  an  ion  at  the  origin 
and/,  which  is  defined  by  Eq.  (16),  is  a  function 
that  reduces  to  the  scattering  amplitude  on  the  en¬ 
ergy  shell  (i.  e. ,  when  k)  =  k%,  =£').  If  we  substi¬ 
tute  Eq.  (16)  in  Eq.  (15)  and  convert  the  sums  to 
integrals,  we  obtain 


TrUio/jGoij  •  ■  ■  G0/n) 

(  1  W  ,3,  f  .3.  eXp[/CEl-  &,!  +  &•  &12+---  +£„•  n)}fi£\  Ci,  E2)  -  -  -  /(£’,  E„,lci) 

°v 2P)r  *'•••] d  - (£*-^-4?"..-^*-^ - - 


where  fty.i,y  =  K/.i  -R,.  In  the  asymptotic  approxi¬ 
mation,  which  amounts  to  evaluating  each  integral 
to  the  lowest  order  in  1/  kR,  the  partial -integration 
method  described  by  Harrison  gives 

TrGo/jGo/j . . .  Wn  =  e"t'sn  ,  (17) 

y.i  h 

where  -  ,  R,.i,y  i ,  s  is  the  sum  of  all  l,,  8,  is  the 
scattering  angle  at  the  yth  site,  andA'z  =  £'.  The 
function  /(£',  8t)  is  now  evaluated  on  the  energy 
shell  and  is,  therefore,  the  scattering  amplitude 
at  energy  £'  and  angle  8r 

From  Eqs.  (5),  (14),  and  (17)  we  find,  for  the 
«-ion  contribution  to  u 


V=- 


Im 


pit  ha  lad 
directions 


fv-n 

Jo  1*1 


fi£,  e<) 


dE 


(18) 

In  arriving  at  this  expression,  we  have  made  use 
of  Harrison’s  observation  that  in  summing  over  all 
ion  positions  in  an  w-ion  array,  each  diagram  con¬ 
necting  the  n  ions  occurs  n  times.  The  sum  over 
paths  and  directions  means  that  in  an  «th  order 
sum  one  should  include  all  n  -vertex  diagrams  for 
which  successive  scatterers  are  on  different  sites, 
and,  for  n  >  2,  that  the  sum  is  to  include  both  di¬ 
rections,  clockwise  and  counterclockwise,  in  which 
the  path  can  be  traversed. 

The  only  difference  between  this  result  and  that 
obtained  by  Lloyd  and  Oglesby  is  that  their  sum  in¬ 
cludes  only  ring  diagrams,  i.e.,  diagrams  in 
which  all  ions  are  distinct.  For  n  2  4,  however, 
it  is  possible  to  draw  proper  n -vertex  diagrams 
(those  for  which  successive  scatterers  are  on  dif¬ 
ferent  sites)  that  involve  fewer  than  n  ions,  and 
for  which  all  ions  are  therefore  not  distinct.  The 


r - ■ — - — — - 

simplest  example  of  this  is  the  fourth-order  dia¬ 
gram  involving  only  two  ions.  In  this  case  there 
are  four  scattering  paths  corresponding  to  propa¬ 
gation  back  and  forth  along  the  line  joining  the 
pair.  Such  a  term  is,  of  course,  of  fourth  order 
and  therefore  asymptotically  negligible  compared 
to  the  second-order  two-vertex  diagram  for  same 
pair.  It  should  not  be  dropped,  however,  in  a  cal¬ 
culation  carried  to  fourth  or  higher  order. 

111.  APPLICATION  TO  TRANSITION  METALS 


For  nontransition  metals,  where  the  scattering 
amplitude  is  a  slowly  varying  function  of  energy, 
one  can  integrate  Eq.  (18)  by  parts  and,  in  the 
asymptotic  approximation,  drop  the  remaining  in¬ 
tegral  which  is  of  order  ( l/kFsf  where  kF  is  the 
Fermi  wave  number.  The  result  is4’5 


V 


tr  = 
n 


Jr 


E  Re 

paths  tad 
direct  teas 


e1*?* 

5 


n 


/a 


£(g ti  e<) 

h 


(19) 


where  we  have  added  the  superscripts  nr  to  indi¬ 
cate  that  Eq.  (19)  is  valid  for  nonresonant  scatter¬ 
ing  only.  The  reason  one  cannot  use  this  result 
for  resonant  scatterers  is  that  the  term  that  was 
dropped  to  obtain  Eq.  (19)  is  an  integral  over  the 
energy  derivative  of  the  scattering  amplitude,  a 
function  with  strong  peaks  near  the  resonance  en¬ 
ergy.  Thus,  for  transition  mete's  one  must  find 
another  way  of  handling  the  ini*  .-al  in  Eq.  (18). 

For  our  model  of  transition -metal  scattering  we 
assume  a  scattering  amplitude  consisting  of  reso¬ 
nant  and  nonresonant  parts8 


/(£,e)  =/“(£,  e)  +  T(£,  6)  , 


where  the  resonant  part  is  of  the  Breit-Wigner 
form 
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/-(£,«- JMCOS.,. 

(20) 

Here  6,  Is  the  so-called  background  phase  shift,’ 
and  Er  and  A  are  the  energy  and  half -width  of  th 
resonance,  respectively.  The  product  of  scatter¬ 
ing  amplitudes  appearing  in  Eq.  (18)  is  therefore 

n  /(£,  e,) 

i- 1 


=/"U>i)/"($.)  ■  •  •  /"(On)  +r(e1)/"(el!)  •  •  •  /"(e„) 
+  •  •  •  +/r(01)/,'(e2)  •  •  •  r(e„) ,  (21) 


slowly  varying  functions  of  energy.  This  allows  us 
to  approximate  the  fw($)  by  their  values  at  the  res¬ 
onance  energy  and  replace  the  wave  number  in  the 
exponential  by  the  approximation 


£=  kf 


E-Zl, 

2 £, 


)■ 


Equation  (20)  then  becomes 


U  =-  g  (-D- /gliir  5L  /’.(cose,) 

nm  »  *  \  *,  /  n;.,/, 


Xlm/“(em<1) .  •  •  /•'(«„) xjs)  ,  (23) 


where  the  £  dependence  of  all  factors  on  the  right- 
hand  side  is  understood.  The  first  term  in  this 
expansion  contains  nonresonant  factors  and  can 
therefore  be  handled  by  the  partial  integration 
method  described  above.  The  last  term,  which 
corresponds  to  pure  resonance  scattering  at  each 
site,  and  the  other  terms,  which  correspond  hybrid 
scattering,  lead  to  integrals  of  the  general  form 

y.  riurih)  •  rc •  •  /"«,) 

hh--l, 

(22) 

We  now  assume  that  the  resonance  energy  lies  well 
above  the  bottom  of  the  conduction  band  (£,  »  0), 
that  the  resonance  region  is  narrow  (A/£r  «  1), 
and,  as  before,  that  the  nonresonant  factors  are 
- 1 


where 

-(irPisp-H'  <“> 

Here  £„(z)  is  the  exponential  integral, 10  y  *A/2kr, 
z,  =  (A  +  iEr)/2kr,  and  z2  =  [A  -  ,(£,  -  Er)}/2kr. 

Explicit  expressions  for  the  two-  and  three -ion 
potentials  are  easily  obtained  by  noting  the  com  - 
binations  of  factors  involved  in  Eq.  (21),  relabel¬ 
ing  indices  as  needed  in  Eq.  (22),  applying  Eq. 
(23),  and  then  adding  the  nonresonant  part.  The 
resulting  potential  for  a  single  pair  (with  /  =  2)  is 


«*(«)=  ^  Re  [/-(£„  it)f+  I *)Xt(s)  -  ~  Xg(s)\  ,  (25) 

where  R  is  the  interionic  distance.  For  a  single  three-ion  path,  the  potential,  summed  over  both  direc¬ 
tions,  is 


vtfuSz,n3)=  ^  Re  ^  ft  WV*  (ze^X^EM,  h,  «,)  -  ~  e*^X,(s) 


*F2(0„e2,03)  +  ^e6,i>X3(s)F3(8u 

where 

Fi(«u  h,  6,)  =  />t(cos«, )/-(*)/-(«,)  +  c.p.  , 
Fi(6u  «j,  03)  =  Pj(cos«l)/>j(cose2)/"r(e3)  +  c.  p.  , 

3 

&2>  ®j) =  jnC  ^ a(cos  fiy)  . 

i*i 

Here  the  1 1  are  interionic  distances  and  the  8,  are 
scattering  angles  for  the  three -ion  diagram,  and 
c.p.  stands  for  cyclic  permutations  of  the  angles 
8U  fi2,  e3.  The  total  energies,  corresponding  to 
V2  and  v3  of  Eq.  (18),  are  obtained  by  summing 
Eqs.  (25)  and  (26)  over  all  two-  and  three-ion 
paths. 


>2.  e3)) .  (26) 

r - —  — — — - - - 

Some  simplification  of  these  results  should  be 
possible  if  s,  the  round-trip  path  length,  is  large 
enough,  for  then  one  could  retain  only  the  leading 
term  in  the  asymptotic  expansions  of  the  exponen¬ 
tial  integrals,  Em(z,s).  However,  in  a  related 
study  involving  the  calculation  of  charge  density 
oscillations  around  resonant  scattering  impurities, 
Mezei  and  Gruner*  found  that  the  use  of  a  similar 
approximation  can  lead  to  serious  errors  in  both 
the  amplitude  and  phase  of  the  oscillations  at  dis¬ 
tances  less  than  a  few  interatomic  spacings.  It  is 
not  difficult  to  see  that  such  errors  might  result 
for  interionic  potentials  as  well.  This  is  because 
In  a  transition  metal  the  energy  difference  Er  -  Er 
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is  of  order  A,  making  the  absolute  value  of  the  ar¬ 
gument  of  the  second  exponential  integral  in  Eq. 

(24)  roughly  [&/Er)krs .  Since  A/£r  «  1  this  is  Like  - 
ly  to  be  a  rather  small  number,  at  least  for  inter¬ 
actions  involving  the  first  few  nearest -neighbor 
shells.  Thus,  for  the  computations  described  be¬ 
low,  we  used  Eqs.  (24)-(26)  without  further  approx¬ 
imation. 


IV.  THIRD-ORDER  CORRECTIONS  IN  TRANSITION 
METALS 

In  Sec.  II  we  found  that  the  configuration-depen¬ 
dent  part  of  the  band-structure  energy  is,  from 
Eqs.  (5)  and  (14), 


IT  « 


£ 


x  ^  Ini  f  Tr(G0<iGo4  •  •  •  G$t^dE  .  (27) 


By  simply  rearranging  the  order  of  summation  we 
can  rewrite  this  expression  in  the  form  of  a  sum 
over  effective  pairwise  potentials.  Thus 

V=  £  V^.Rj)  ,  (28) 

where 


■U(Kl,R8)  = 

with 


2 

it 


Trtoo^Go^Jdfi  , 


(29) 


it-3  n 


£  ' "  ’  £  G0faG()<4  •  •  •  G0t„ 


(30) 


Except  in  the  lowest  order  approximation,  fg»fg, 
is,  of  course,  not  a  true  pairwise  poten¬ 
tial  in  the  usual  sense,  because  it  includes  inter¬ 
actions  with  all  other  ions  in  the  crystal.  On  the 
other  hand,  the  effective  potential  formalism  is 
useful  for  our  present  purpose  because  it  shows 
how  one  should  add  multi-ion  corrections  to  the 
isolated-pair  approximation  to  obtain  the  correct 
result  for  the  band -structure  energy  written  as  a 
sum  of  pairwise  interactions.  Thus,  for  example, 
we  find  that  in  the  asymptotic  approximation  to 
third  order 


U®j,  5g)  -  U4(  |  Rg  -  | )  +  t  £  Uj®j,  ftj) 

5Sfi 


(31) 


where  vt  and  Vs  are  given  by  Eqs.  (25)  and  (28). 
We  note  once  again  that  the  correction  term  is  de¬ 
pendent  on  the  positions  of  ions  1  and  2  with  re¬ 
spect  to  the  other  ions  and  is  not,  therefore,  a 


simple  function  of  separation  distance. 

For  our  initial  numerical  investigation  of  the 
role  of  multi -ion  interactions  in  transition  metals 
we  applied  Eq.  (31)  to  the  calculation  of  nearest- 
neighbor  interactions  in  a  face -centered -cubic 
crystal.  The  nearest -neighbor  distance  was  taken 
to  be  5  a.  u. ,  the  resonance  parameters  were  those 
used  by  Pettifor, 11  namely,  £r  =  0.5  Ry  and  A  =  0.03 
Ry,  and  pure  s-wave  scattering  with  a  phase  shift 
of  0.1  rad  was  assumed  for  the  nonresonant  part. 
To  simulate,  in  a  rough  way,  the  behavior  of  the 
effective  potential  in  a  transition -metal  series,  we 
allowed  the  Fermi  energy  to  vary  from  0.4  to  0.6 
Ry.  The  sum  over  neighboring  ions  included  only 
the  nearest  neighbors  of  ions  1  and  2.  While  this 
cluster  is  too  small  for  accurate  quantitative  pre¬ 
dictions,  density -of -states  calculations  by  a  simi¬ 
lar  cluster  method14, ls  indicate  that  the  nearest- 
neighbor  model  is  adequate  for  qualitative  studies 
of  the  effects  of  the  splitting  of  the  d  resonance  in¬ 
to  subbands.  Thus,  our  calculations  should  provide 
a  reliable  indication  of  at  least  one  aspect  of  the  ef¬ 
fects  of  multi-ion  interactions,  i.e.,  the  role  of  d~ 
band  splitting  in  the  effective  interionic  interac¬ 
tion. 

The  results  are  plotted  in  Fig.  1  as  a  function 
of  the  position  of  the  Fermi  level  with  respect  to 
the  resonance  energy.  The  main  features  in  the 
isolated-pair  potential,  shown  as  a  dashed  curve  in 
this  figure,  are  the  strong  attractive  peak  that  oc¬ 
curs  as  the  Fermi  level  passes  through  the  reso- 


FIG.  1.  Band-structure  contributions  to  the  interionic 
potential  for  nearest  neighbors  in  a  hypothetical  face- 
centered-cubic  transition-metal  series.  The  dashed 
curve  is  the  potential  for  an  isolated  pair,  plotted  as  a 
function  of  the  position  of  the  Fermi  level  with  respect 
to  the  resonance  energy,  while  the  solid  curve  was  ob¬ 
tained  by  adding  three-ion  scattering  corrections  as 
described  in  the  text.  These  results  indicate  that  multi- 
ion  interactions  play  an  essential  role  in  determining 
the  effects  of  rf-band  splitting  on  interionic  interactions 
in  transition-metal  crystals. 
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nance  and  the  weaker  repulsive  peak  at  higher  en¬ 
ergies.  However,  the  solid  curve,  which  was  ob¬ 
tained  by  adding  third-order  corrections,  shows 
that  the  isolated-pair  behavior  is  significantly  al¬ 
tered  by  multi -ion  interactions  in  a  solid.  It  would 
appear  from  these  results  that  multi-ion  scattering 
tends  to  shift  the  attractive  peak  to  a  lower  energy 
while  introducing  a  new  repulsive  component  at  an 
energy  slightly  above  resonance.  This  behavior 
corresponds,  at  least  in  a  qualitative  way,  to  the 
development  of  an  additional  peak  in  a  cluster  cal¬ 
culation  of  the  d-band  density  of  states. 12,11  Ac¬ 
cording  to  the  localized  orbital  model, 14  these 
peaks  can  be  associated  with  bonding  states  at  low¬ 
er  energies  and  antibonding  states  at  higher  ener¬ 
gies.  The  predicted  nature  of  the  peak  structure 
of  the  effective  interionic  potential,  i.e.,  attrac¬ 
tive  at  lower  energies  and  repulsive  at  higher  en¬ 
ergies,  is,  of  course,  in  agreement  with  the  bond - 
ing-antibonding  orbital  picture. 

While  these  results  are  encouraging,  because 
they  show  that  a  relatively  simple  third -order  cor¬ 
rection  gives  an  interionic  potential  in  qualitative 
accord  with  what  is  known  about  the  electronic 
structure  of  transition  metals,  it  must  be  remem¬ 
bered  that  a  number  of  approximations  are  in¬ 
volved.  The  very  fact  that  third-order  terms  are 
important,  for  example,  makes  one  wonder  about 
the  effects  of  fourth-  and  higher-order  terms.  The 
truncation  of  the  cluster  expansion  at  nearest- 
neighbor  shells  is  also  a  matter  of  concern,  par¬ 
ticularly  since  such  small  cluster  approximations 
are  known  to  lead  to  poor  results  for  the  free-elec- 
tron  and  hybridization  contributions  to  the  density 
of  states. 11 

There  are,  in  addition,  some  unanswered  ques¬ 
tions  at  a  more  fundamental  level.  The  validity  of 
the  independent-electron  self-consistent -potential 
model,  which  is  implicit  throughout  the  develop¬ 
ment,  is  itself  subject  to  some  criticism  in  the 
present  context.  For  example,  it  is  not  clear  that 
one  can  in  fact  find  an  approximate  energy -inde¬ 
pendent  exchange  potential  [such  as  is  required  by 
Eq.  (8)]  that  will  lead  to  an  adequate  description  of 
nearest-neighbor  interactions  near  the  noble  metal 


end  of  a  transition-metal  series.11  [In  a  multi-ion 
treatment  of  the  effective  pairwise  potential,  near¬ 
est-neighbor  interactions  are  involved,  even  in  an 
asymptotic  theory;  see,  for  example,  Eq.  (31).) 
Furthermore,  even  if  we  ignore  this  matter  by  as¬ 
suming  that  a  suitable  exchange  approximation  can 
be  found,  there  remains  the  problem  of  removing 
half  the  electron-electron  interaction  energy, 
which  is  always  counted  twice  in  a  calculation  of 
the  total  energy  in  the  independent -electron  approx¬ 
imation.  While  it  is  rather  easy  to  subtract  this 
extra  energy  in  a  second -order  pseudopotential 
treatment, 1 2 3 4 5 * 7  we  know  of  no  practical  way  to  proceed 
within  the  framework  of  the  /-matrix  formalism. 
Thus,  when  treating  transition  metals,  where  res¬ 
onance  scattering  demands  a  /  -matrix  (or  equiva¬ 
lent)  description  of  scattering  by  individual  ions, 
it  is  common  practice  to  simply  ignore  the  double - 
counting  of  electron-electron  energies, 2,5  as  we  did 
here.  In  the  final  analysis,  what  all  of  this  really 
means  is  that  in  order  to  develop  a  tractable  mod¬ 
el,  we  have  assumed  that  the  approximations  dis¬ 
cussed  above  are  of  secondary  importance,  and 
that  the  major  features  of  multi -ion  interactions  in 
transition  metals  are  adequately  accounted  for  by 
a  simple  multiple  scattering  calculation  for  a 
small  cluster  of  ions . 

There  is,  therefore,  ample  reason  to  question 
the  quantitative  accuracy  of  such  features  as  the 
positions,  heights  and  widths  of  the  peaks  shown 
in  Fig.  1.  On  the  other  hand,  the  results  of  other 
calculations  and  comparisons  with  experiment11-11 
indicate  that  our  calculational  model  should  provide 
a  reasonable  description  of  resonant  scattering  ef¬ 
fects,  and  that  one  should  expect  the  main  result, 
the  existence  of  attractive  and  repulsive  structures 
in  the  effective  interionic  potential,  to  persist  in 
more  accurate  calculations. 
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The  theory  of  multt-ion  interactions  is  applied  to  calculations  of  the  stacking-fault  energies  of  transition 
metals  It  is  found  that  observed  trends  in  the  stabilities  of  close-packed  transition  metals  can  be  explained  on 
the  basis  of  a  simple  model  in  which  resonant  scattering,  three-ion  interactions  dominate. 


I  INTRODUCTION 

The  results  of  recent  calculations  of  stacking- 
fault  energies  and  comparisons  with  experimental 
data  indicate  that  a  pseudopotential  expansion 
carried  to  second  order  in  the  energy  is  adequate 
for  the  treatment  of  nontransition  metals. 1,2  As 
Heine  and  Weaire3  had  noted  earlier,  it  was  to  be 
expected  that  a  second-order  theory  would  be  suc¬ 
cessful  in  this  particular  application  because  the 
pseudopotential  matrix  elements  involved  in  the 
perturbation  expansion  are  small. 

For  transition  metals,  on  the  other  hand,  it  is 
unlikely  that  such  a  simple  theory  will  suffice  be¬ 
cause  the  resonant  scattering  of  d  electrons  by 
transition- metal  ions  leads  to  large  matrix  ele¬ 
ments  and,  hence,  to  slower  convergence  of  the 
perturbation  expansion.  Indeed,  even  for  the 
noble  metals,  where  one  might  hope  that  the  ef¬ 
fects  of  resonant  scattering  are  not  too  pronounced, 
calculated  stacking- fault  energies  based  on  a  sec¬ 
ond-order  approximation  are  in  poor  agreement 
with  experimental  results. 1,2  It  seems  clear, 
therefore,  that  a  successful  theory  of  stacking- 
fault  energies  in  transition  metals,  and  perhaps 
the  noble  metals  as  well,  must  include  the  effects 
of  higher-order  terms  in  the  perturbation  expan¬ 
sion.  4  The  purpose  of  the  work  reported  here 
was  to  see  if  an  extension  to  third  order  will  ade¬ 
quately  account  for  observed  trends  in  the  stack- 
ing-fault  energies  and  relative  stabilities  of  close- 
packed  transition  metals. 

In  analogy  with  Harrison’s  formulation  of  the 
second-order  theory5  of  stacking- fault  energies, 
one  might  attempt  a  third-order  calculation  by 
making  use  of  the  formal  expression*  for  the  third- 
order  energy  in  terms  of  the  plane-wave  matrix 
elements  of  the  pseudopotential  and  the  structure 
factors  for  the  perfect  and  faulted  crystal.  It  is 
well  known,  however,  that  calculations  of  third- 
order  energies  by  this  method  are  extremely 
complex,  even  for  perfect  nontransition-metal 
crystals.  7,8  Since  the  prospect  of  extending  such 
a  calculation  to  stacking  faults  with  the  additional 
complication  of  resonant  scattering  is  indeed  for¬ 
midable,  for  the  investigation  reported  here  we 

13 


choose  an  alternative,  approximate  approach 
based  on  the  theory  of  multi-ion  interactions. 9,10 

In  short,  our  approach  is  based  on  Harrison’s 
observation  that  the  third-order  term  in  the  per¬ 
turbation  expansion  of  the  total  energy  can  be 
written  as  a  sum  of  three- ion  interaction  energies. 9 
The  third-order  term  in  the  expression  for  the 
stacking-fault  energy  is  therefore  the  difference 
between  the  three-ion  sum  for  the  faulted  crystal 
and  the  corresponding  sum  for  the  perfect  crystal. 
By  truncating  these  sums,  i.  e. ,  by  keeping  only 
those  three-  ion  interactions  judged  to  be  dominant 
in  the  determination  of  the  stacking-fault  energy, 
we  obtain  an  approximation  to  the  third-order  en¬ 
ergy. 

In  the  text  of  this  paper  we  will  first  describe 
in  more  detail  the  method  for  calculating  the 
third-order  energy.  An  approximate  method  for 
estimating  the  second-order  energy  will  then  be 
derived,  followed  by  a  discussion  of  the  electron- 
ion  scattering  model  used  in  numerical  computa¬ 
tions.  The  results  of  our  study,  which  are  pre¬ 
sented  in  the  concluding  section,  show  that  for  most 
transition  metals  it  is  the  three-ion  terms,  not 
pairwise  interactions,  that  comprise  the  dominant 
contribution  to  the  stacking-fault  energy.  Al¬ 
though  quantitative  agreement  with  experimental 
data  on  stacking-fault  energies  is  only  fair,  the 
calculations  do  adequately  explain  observed  trends 
in  the  relative  stabilities  of  close-packed  struc¬ 
tures  for  the  first  three  transition- metal  series. 

II.  CALCULATION  OF  THE  THREE  ION  ENERCY 

Our  starting  point  is  the  following  general  ex¬ 
pression  for  the  three- ton  contribution  to  the 
stacking-fault  energy: 

r^EEE 

(i) 

where  Rf  and  Rf  are  ion  positions  in  the  faulted 
(F)  and  perfect  (P)  crystals,  A  is  the  fault  area, 
and  Uj  is  the  three-ion  interaction  energy  defined 
in  Ref.  10.  Since  t)s  already  includes  a  sum  over 
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cyclic  permutations  of  ion  positions  (i.  e. ,  a  sum  of 
identical  three-ion  diagrams)  only  distinct  diagrams 
are  to  be  included  in  the  sums  on  i,  j,  and  k.  Also, 
terms  in  which  all  three  ion  positions  lie  on  the 
same  side  of  the  fault  plane  do  not  contribute  be¬ 
cause  in  such  cases  the  three-ion  energies  are  the 
same  in  the  faulted  and  perfect  crystals.  Thus  the 
only  diagrams  that  need  be  considered  are  those 
distinct  diagrams  where  two  of  the  position  vectors 
terminate  on  opposite  sides  of  the  fault  plane.  We 
therefore  choose  R,  to  be  the  position  of  an  ion  on 
one  side  of  the  fault  plane  and  R,  a  position  on  the 
opposite  side,  with  R,  on  either  side. 

Figure  1  is  an  example  of  a  pair  of  three  ion  con¬ 
figurations  that  give  a  nonvanishing  contribution  in 
the  calculation  of  the  intrinsic  fault  energy  in  an  fee 
crystal.  Here  R,  and  Ry  are  the  same  in  the  per¬ 
fect  and  faulted  structures  while  R^  terminates  at 
the  positions  indicated  by  Rf  and  R'  in  the  perfect 
and  faulted  crystals,  respectively. 

For  this  particular  geometrical  arrangement,  it 
is  easily  seen  that  for  each  ion  position  Ry  below 
the  fault  plane  there  are  three  equivalent  diagrams 
corresponding  to  the  three  nearest -neighbor  posi¬ 
tions  of  R,.  Also,  there  are  three  more  equivalent 
diagrams  in  the  mirror-image  configuration  in 
which  R,  and  Ry  terminate  above  the  plane  while  R, 
lies  below.  Thus,  assuming  that  the  energy  dif¬ 
ference  in  Eq.  (1)  is  the  same  for  all  Ry  below  the 
fault  plane,  we  obtain,  for  that  part  of  the  stacking 
fault  energy  due  to  the  interactions  illustrated  in 
Fig.  1, 

Ay,=( l/2w)(6vfBA-6VfBC)  , 

where  io  is  the  area  per  ion  and  and  V*”0  are 
the  energies  corresponding  to  the  particular  dia¬ 
grams  considered  here.  Similar  arguments  can, 


FIG.  1.  Typical  pair  of  three-ton  interactions  that  con¬ 
tribute  to  the  intrinsic  stacldng-fault  energy  for  face- 
centered -cubic  crystals.  Here  SJ  and  8J  are  ion  posi¬ 
tions  in  the  perfect  (P)  and  faulted  (F)  crystals,  while  R,  and 
Ry  are  the  same  in  both  caseB. 


foult  plan* 


FIG.  2.  Schematic  illustration  of  three-ion  interac¬ 
tions.  Terms  included  in  the  truncated  stacldng-fault 
energy  sum  were  those  for  which  either  a  $  30°  or  R 
4 20  a.u. 


of  course,  be  applied  to  other  three-ion  configura¬ 
tions  thus  leading  to  the  formula 

ya  =  -^;(E  tffvf (fl  -  £  ATflfU))  ,  (2) 

where  N*  is  the  number  of  equivalent  diagrams  of 
the  fth  geometry  in  the  faulted  crystal.  v'(i)  is 
the  corresponding  energy,  the  IV f  andT'f(f)  are  the 
diagram  weights  and  energies  for  the  perfect  fee 
structure,  and  the  sum  is  over  all  nonequivalent 
three-ion  diagrams. 

In  the  applications  discussed  in  Sec.  V,  the  only 
diagrams  considered  were  those  for  which  at  least 
two  of  the  three  ions  were  nearest  neighbors  and 
for  which  the  third  ion  was  in  either  the  adjacent 
Or  next-nearest  plane,  as  illustrated  in  Fig.  2. 

In  addition,  the  sums  over  diagrams  were  limited 
to  terms  for  which  either  a«30°  orfls  20  a.u. 

(see  Fig.  2).  While  the  truncation  of  the  sums  at 
these  particular  values  of  a  and  R  is  somewhat 
arbitrary,  it  can  be  seen  from  the  expression  for 
V,  given  in  Ref.  10  that  the  three-ion  energy  is  in¬ 
versely  proportional  to  the  product  of  the  three  in¬ 
terior  path  lengths  times  the  roundtrip  distance, 
thus  causing  Vj  to  fall  off  rather  rapidly  with  in¬ 
creasing  R.  Also,  because  the  resonant  scattering 
terms  in  t'j  contain  products  of  the  Legendre  poly¬ 
nomials  P,  (cos©,),  where  ©,  is  the  scattering  an¬ 
gle  at  the  tth  ion  site,  the  energy  tends  to  be  larg¬ 
er  for  diagrams  with  a*  0  than  for  diagrams  with 
larger  a.  Actual  computations  of  V,  confirmed 
these  expected  trends,  in  that  they  showed  that  en¬ 
ergies  corresponding  to  the  larger  values  of  a  and 
R  are  small,  usually  about  10%  of  the  dominant 
terms.  Still,  it  should  be  noted  that  in  some  cases 
the  larger  terms  tend  to  cancel  and,  when  this  hap¬ 
pens,  the  truncation  error  can  be  significant.  We 
will  return  to  this  point  when  discussing  the  results 
of  our  computations  in  Sec.  V. 

III.  CALCULATION  OF  THE  PAIRWISE  INTERACTION 
ENERGY 

To  calculate  that  part  of  the  stacking  fault  en¬ 
ergy  due  to  pairwise  Interactions  we  used  the  for- 
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malism  of  Blandin  et  al. 11  They  showed  that  the 
pairwise  contribution  to  y  can  be  written 

■O 

y2  =  X!  M»)A0(«A)  . 

n-2 

where  h  is  the  distance  between  close-packed 
planes,  M«)  is  an  integral  weight  corresponding 
to  a  particular  fault  configuration,  and  A <f>(nh)  is 
an  interplanar  potential  difference. 1,1  Blandin 
et  al. ,  also  showed  that  the  potential  difference 
is  related  to  the  energy-wave  number  character¬ 
istic  4>(c/)  as  follows: 

nk)*irJ. Z  +*2>1/2i  .  <3> 

where  g  is  the  magnitude  of  the  smallest  nonvan¬ 
ishing  reciprocal-lattice  vector. 

From  the  Fourier  transform  relationship1  be¬ 
tween  $(<?)  and  the  pairwise  interaction  energy  u2, 
we  find,  using  the  asymptotic  approximation, 10 

VZ(X)~-—  fEF Im  dE 

IT  •'q  ti 

that  for  q>  2 kf,  where  kt  is  the  Fermi  wave  num¬ 
ber, 

*{q)*  ■  wa  C 1/{E’ 12)1 2  (*  cos2Hi^l 

+  |  sin2rj^  dE  (4) 

where  fi0  is  the  volume  per  ion,  f(E,  n)  is  the  scat¬ 
tering  amplitude  at  energy  E  and  scattering  angle 
ir,  and  k  =  V"E .  The  angle  tj  is  defined  in  terms  of 
the  scattering  phase  shifts  5,  by  tan?j  =  a/0,  where 

a  =  £(- l)'(2f+ l)cos6,sin6,  , 

i 

0=]C(_  1)'(2/+  l)sin*8,  . 

i 

The  condition  9>2x,-is  satisfied  here  because, 
from  Eq.  (3),  q>  g  and  g>lKr  for  all  transition 
metals.  Substitution  of  Eq.  (4)  in  Eq.  (3)  gives,  to 
first  order  in  l/tih  and  neglecting  terms  of  order 
exp  (~g>ih) 

\f(E,n) |2cos2r?-—^ — dE  , 

where 

»,  =  l  Ugf-E)"1. 

In  applying  this  result  in  the  computations  described 
below  the  final  integration  was  numerically  evalu¬ 
ated,  and  only  the  term  corresponding  to  next- 


nearest-neighbor  plane  interactions  (»  =  2)  was  re¬ 
tained  in  the  stacking  fault  energy  sum. 

IV.  CALCULATION  OF  PHASE  SHIFTS 

Our  choice  of  a  model  for  the  calculation  of  scat¬ 
tering  phase  shifts  for  transition  metal  ions  is 
based  largely  on  the  work  of  Pettifor. 12  He  showed 
that  one  can  reproduce,  with  reasonable  accuracy, 
the  results  of  more  detailed  calculations  of  the  den¬ 
sities  of  states  of  transition  metals,  by  basing  a 
simpler  calculation  on  the  assumption  that  the  d- 
electron  resonance  energy  Er  and  width  A  are  con¬ 
stants  for  the  first  three  transition-metal  series. 

In  the  calculation  reported  here  we  made  the  ad¬ 
ditional  approximations  that  the  nonresonant  scat¬ 
tering  phase  shifts  can  be  derived  from  a  pseudo¬ 
potential  which,  again,  is  the  same  for  all  ele¬ 
ments,  and  that  the  nearest-neighbor  distance  is 
the  same  (5.0  a.u. )  for  all  elements.  This  leaves 
us  with  a  rather  simple  calculational  model  in 
which  the  only  parameter  that  distinguishes  one 
element  from  the  next  is  the  Fermi  energy.  How¬ 
ever,  except  for  the  approximation  concerning 
phase  shifts  for  nonresonant  scattering,  this  is  the 
same  model  used  by  Pettifor  in  his  calculations  of 
the  relative  energies  of  the  hep,  fee,  and  bcc 
structures.  The  fact  that  Pettifor’s  results  are  in 
accord  with  the  observed  structures  of  transition 
metals  suggests  that  the  model,  through  obviously  an 
idealization,  forms  a  reasonable  basis  for  the  study 
of  stacking  fault  energies  as  well. 

For  the  resonance  parameters  we  used  Moriar- 
ty’s  values  for  copper13  (Er  =  0.  33  Ry  and 
A  =  0.014Ry).  Fermi  energies  as  a  function  of 
valence  Z  were  determined  by  adding  an  average 
of  Pettifor’s  calculated  values  of  Ef~Er  for  the 
fee  and  hep  structures  to  the  value  chosen  for  Er. 

The  value  of  the  resonant  part  of  the  d  electron 
phase  shift  was  determined  from  the  formula 

tan(8,-6*)  =  A/(£r-E)  , 

where  62  and  6t  are  the  resonant  and  nonresonant 
parts  of  the  l  =2  phase  shift 

To  calculate  the  nonresonant  scattering  phase 
shifts  we  used  the  empty  core  potential,  the  core 
radius  being  that  given  by  Ashcroft  and  Langreth1’ 
for  copper.  The  phase  shifts  were  calculated  in 
the  first  Born  approximation  and  were  based  on  the 
uniform  screening  charge  assumption. 15  Because 
there  is  considerable  uncertainty  as  to  the  validity 
of  such  a  simple  model,  we  performed  two  sets  of 
calculations  using  effective  valences  {Zs)  of  one  and 
two  in  the  nonresonant  scattering  pseudopotential. 

V.  RESULTS  AND  DISCUSSION 

The  results  of  our  calculations  are  shown  in  Fig. 

3.  These  plots  show  the  stacking  fault  energy  for 
intrinsic  faults  in  fee  crystals  as  a  function  of  Z , 
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VALENCE 


FIG.  3.  Calculated  intrinsic  stacking-fault  energies  for 
face-centered-cubic  crystals.  The  parameter  Zs  is  the 
valence  assumed  in  the  calculation  of  nonresonant  scat¬ 
tering  phase  shifts.  The  dashed  parts  of  the  curves 
correspond  to  values  of  the  total  valence  (Including  d 
electrons)  where  neither  close-packed  structure  is  sta¬ 
ble.  Experimental  data,  shown  here  by  the  symbol  x , 
are  taken  from  Ref.  16  for  copper  (Z  =  11)  and  nickel-co¬ 
balt  alloys  (9  <  Z  <  10),  and  from  Ref.  17  for  cobalt 
<Z  =  9), 


the  total  valence,  for  both  values  of  Z,.  The 
dashed  part  of  the  curve  corresponds  to  those  val¬ 
ues  of  Z  where  neither  close-packed  structure  is 
stable  and  where,  therefore,  no  comparisons  with 
observed  stable  structures  or  stacking  fault  ener¬ 
gies  are  possible.  Also  shown  are  experimental 
data  for  cobalt  (Z  =  9),  nickel  (Z  =  10),  copper 

=11),  and  cobalt-nickel  alloys.  The  stacking 
fault  energy  for  hep  cobalt  is  shown  as  a  negative 
number  because  the  theoretical  expression  for  the 
energy  of  an  intrinsic  fault  in  an  hep  crystal  is  ap¬ 
proximately  equal  to  minus  the  expression  for  the 
intrinsic  fee  energy.  As  can  be  seen  from  the 
stacking  sequences  for  the  perfect  and  faulted  crys¬ 
tals,  this  relationship  is  exact  if  one  ignores  all 
interactions  involving  third  and  more  distant  neigh¬ 
bors  planes. 

The  first  point  to  be  noted  here  is  that  the  two 
calculated  curves  are  in  reasonably  good  agree¬ 
ment  for  most  values  of  Z.  Since  these  two  curves 
are  based  on  nonresonant  scattering  phase  shifts 
that  differ  by  about  a  factor  of  2,  this  must  mean 
that  the  details  of  nonresonant  scattering  are  rela¬ 
tively  unimportant,  and  that  the  general  trend  of 
stacking  fault  energy  versus  valence  is  determined 
largely  by  resonant  scattering  properties.  It 
should  be  noted,  however,  that  the  resonance  width 
A,  which  we  have  taken  from  Moriarty’s  calcula¬ 


tion,  is  related  to  nonresonant  scattering  proper¬ 
ties  through  hybridization.  One  would  therefore 
expect  that  a  first-principles  calculation  would 
show  a  stronger  dependence  on  nonresonant  scat¬ 
tering  properties  than  is  indicated  here. 

Another  point  to  be  noted  regarding  Fig.  3  is  that 
the  stability  of  the  fee  phase  against  faulting  is  cor¬ 
rectly  predicted  for  most  values  of  Z.  Thus,  for 
2>9  we  obtain  positive  values,  which  indicates 
stability  against  faulting  in  fee  crystals,  while  for 
6  <  Z  <  9  and  Z  <  4,  where  the  stable  close-packed 
structure  is  hep,  we  obtain  negative  fault  energies. 
The  only  notable  exception  is  the  Z,  =  2  curve  at 
Z  =  4  (titanium,  zirconium,  and  hafnium),  where 
the  calculation  indicates  that  the  fee  structure  is 
stable.  This  failure  may  well  be  due  to  the  approx¬ 
imate  treatment  of  nonresonant  scattering,  since 
the  energy  for  Zs  =  1  has  the  correct  sign.  On  the 
other  hand,  the  fact  that  the  absolute  value  stacking 
fault  energy  is  much  smaller  at  Z  =  4  than  at  other 
values  of  Z  indicates  that  there  is  considerable  can¬ 
cellation  in  the  sum  of  two-  and  three-ion  interac¬ 
tion  energies,  and  that  the  truncation  error  men¬ 
tioned  above  may  therefore  be  significant.  It 
should  also  be  noted  that  a  similar  situation  exists 
at  Z  =  9  (cobalt ,  rhodium,  and  iridium).  Here, 
however,  one  of  the  elements  (cobalt)  is  in  fact 
stable  in  the  hep  phase,  while  the  other  two  form 
stable  fee  crystals. 

It  is  of  interest  to  compare  these  results  with 
Pettifor’s  calculations  of  the  relative  energies  of 
close-packed  phases.1*  As  was  noted  previously, 
our  assumptions  of  constant  resonance  parameters 
and  nearest-neighbor  distance  are  consistent  with 
Pettifor’s  model  although  he  used  different  numer¬ 
ical  values  for  the  resonance  parameters.  How¬ 
ever,  the  principal  difference  between  his  model 
and  ours  lies  in  the  treatment  of  nonresonant  scat¬ 
tering.  In  spite  of  these  differences  our  results 
for  the  relative  stabilities  of  the  fee  and  hep  struc¬ 
tures  as  a  function  of  Z  are  in  reasonably  good 
agreement  with  Pettifor’s.  Thus  he  finds  that  there 
is  a  region  of  hep  stability  near  Z  =  4  followed  by  an 
fee  stable  region  near  Z  =  5,  an  fee  to  hep  transi¬ 
tion  near  Z-  7  and  finally,  an  hep  to  fee  transition 
near  Z  =  9.  We  take  this,  and  the  fact  that  the  re¬ 
sults  agree  with  observed  trends  in  crystal  struc¬ 
ture,  as  indications  that  our  simple  scattering 
model  is  adequate  as  a  first  approximation,  and 
that  we  have  included  the  most  significant  terms  in 
our  approximate  summation  of  the  multi-ion  ex¬ 
pansion. 

Regarding  the  calculated  values  of  the  stacking- 
fault  energies,  as  can  be  seen  in  Fig.  3,  agree¬ 
ment  with  experimental  data  is  satisfactory  for 
copper  and  cobalt,  but  poor  for  nickel  and  cobalt - 
nickel  alloys.  In  view  of  the  very  simple  model 
that  was  used  in  calculating  electron-ion  scattering 
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FIG.-l.  Pairwise  (>2)  and  three-ion  (y3)  contributions  to 
the  stacking  fault  energy  for  2,  =  1.  At  almost  all  values 
of  the  valence  the  three-ion  energy  is  much  greater  than 
the  pairwise  contribution,  which  is  multiplied  by  10  in 
this  plot. 


phase  shifts,  we  do  not  consider  these  comparisons 
particularly  significant,  except,  perhaps,  as  an 
indication  that  a  more  careful  treatment  is  needed 
for  quantitative  comparisons  with  experiment. 

Finally,  in  Fig.  4  we  show  the  two-  and  three- 
ion  contributions  to  the  stacking  fault  energy  for 
Z ,  I  (the  results  for  Z,~  2  are  similar).  The 
important  point  to  be  noted  here  is  that  three-ion 
interactions  dominate  at  almost  all  values  of  Z. 
Although  one  might  expect,  on  the  basis  of  our  ear¬ 
lier  calculation  for  the  noble  metals, 1  that  a  more- 
accurate  nonasymptotic  calculation  would  yield 
larger  values  for  the  pairwise  contribution,  we 


still  believe  that  the  dominance  of  three-ion  terms 
exhibited  here  is  at  least  qualitatively  correct. 

One  reason  for  this  is  that  the  pairwise  energies 
are  so  small  compared  to  three-ion  energies  that 
the  errors  introduced  through  the  asymptotic  ap¬ 
proximation,  being  roughly  of  the  same  order  of 
magnitude  as  the  pairwise  predictions  themselves, 1 
are  expected  to  be  insignificant  compared  to  the 
three-ion  energies. 

Another  reason  is  that  from  Peitifor’s  calcula¬ 
tion,  and  from  the  observed  stability  of  close- 
packed  phases  as  a  function  of  Z ,  one  would  expect 
three  hep- fee  transitions  in  a  transition  metal  ser¬ 
ies.  The  three-ion  energy  does,  in  fact,  show 
three  such  transitions  while  the  pairwise  energy 
has  only  two.  Thus  assuming  only  that  the  shapes 
of  the  two-  and  three-ion  curves  shown  in  Fig.  4 
are  correct,  one  would  expect  three-ion  contribu¬ 
tions  to  dominate. 

In  conclusion,  therefore,  we  have  demonstrated 
that  observed  trends  in  the  stabilities  of  close- 
packed  transition  metals  can  be  explained  on  the 
basis  of  a  third-order  calculation  involving  a  sim¬ 
ple  resonant  scattering  model  of  electron-ion  in¬ 
teractions.  The  calculations  indicate  that  such 
trends  are  governed  largely  by  three-ion  interac¬ 
tions  at  most  values  of  the  valence.  Comparisons 
with  experimentally  measured  stacking- fault  ener¬ 
gies  show,  as  expected,  that  a  more  careful  treat¬ 
ment  is  needed  for  quantitative  predictions  of  stack¬ 
ing-fault  energies.  Still,  the  degree  of  success 
realized  in  the  prediction  of  structural  trends  and, 
at  least,  the  correct  sign  and  order  of  magnitude 
of  stacking-fault  energies,  supports  our  principal 
conclusions  regarding  the  dominant  roles  of  reso¬ 
nant  scattering  and  three-ion  interactions. 
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It  is  pointed  out  that  the  Lloyd-Sholl  expression  for  the  total  energy  is  self  consistent  through  third  order  in  the 
electron -ion  pseudopotential. 


Several  years  ago  Lloyd  and  Shod  [1  ]  and  later 
Brovman  et  al.  [2]  independently  developed  an  ex¬ 
pression  for  the  third  order  term  in  the  expansion  of 
the  total  energy  of  a  metal.  Although  self-consistency 
of  the  electron -ion  pseudopotential  is  not  discussed 
in  these  works,  proper  consideration  of  nonlinear 
screening  to  second  order  in  the  pseudopotential, 
which  is  required  in  a  self-consistent  third  order  theory 
of  the  total  energy,  is  implicit  in  the  calculations.  An 
explicit  expression  for  the  second  order,  self-consis¬ 
tent  pseudopotential  was  presented  later  by  Bertoni 
et  al.  (3 1 ,  but  they  did  not  show  how  the  second 
order  screening  term  enters  the  total  energy  expres¬ 
sion.  The  purpose  of  this  letter  is  to  clarify  the  role  of 
nonlinear  screening  in  the  third  order  energy.  We  show 
that  the  third  order  terms  associated  with  nonlinear 
screening  cancel  in  the  total  energy  calculation  and 
that  the  Lloyd-Sholl  result  is  therefore  self-consistent 
to  third  order. 

The  total  energy  of  a  crystal  can  be  written  as  the 
sum  of  a  free  electron  energy,  a  direct  interionic  inter¬ 
action  energy,  and  the  band  structure  energy,  which 
is  the  structure  dependent  part  of  the  electron -ion 
interaction  energy.  Because  screening  effects  enter 
only  in  the  band  structure  energy,  we  need  not  con¬ 
sider  the  other  terms.  According  to  Lloyd  and  Shell 
the  band  structure  energy  per  ion  is 


y.  1  -e(P) 
2  p  *  o  V(p) 


\Wp\2 


£o  £ 

2  p  *  o  n0m 


+  n0  S  ,P2'Pi)wp,  wPl  ftp, 

p\pi  (i) 

+  higher  order  terms, 

where  Wp  is  the  self-consistent  pseudopotential,  efp) 
is  the  dielectric  function,  and  is  the  Coulomb 

interaction  corrected  for  exchange  and  correlation 
energies.  It  is  related  to  the  corrected  dielectric  func¬ 
tion  e(p)  and  the  Hartree  dielectric  function  e0(p),  as 
follows: 


n0^)  =  ^ 

p2 


1  -  e(P) 

•  -  eo(P) ' 


The  potential  is  that  due  to  valence  electrons  and 
g(pi ,p2,Pj)is  the  function  defined  by  Lloyd  and 
Sholl,  with  py  =  -  (/>j  +  p2).  The  second  term  in  eq. 
(1)  is  the  correction  for  double  counting  of  electron 
-electron  interactions  in  a  self-consistent  field  treat¬ 
ment. 

In  the  linear  screening  approximation  the  sel  Con¬ 
sistent  pseudopotential  is  given  by 

%  ~  W”)  =  W™le{p), 

where  Hjj°n  is  the  ionic  pseudopotential.  To  account 
for  higher  order,  nonlinear  screening,  we  take 


Wp*W<"+AH>p, 


(2) 
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where  AW  is  a  correction  for  second  and  higher  order 
terms  in  the  screening  charge  density  expansion.  The 
valence  electron  potential  is  therefore 

w;  =  W?  >  -  W™  +  A  Wp  -  +  A  V 

If  we  substitute  eqs.  (2)  and  (3)  in  the  first  two 
terms  in  eq.  (1)  we  find  that  the  terms  linear  in  AWp 
cancel  exactly.  This  leaves  only  a  term  in  I AWp\2, 
which  is  of  fourth  order,  plus  higher  order  corrections 
from  the  other  terms  in  eq.  (1).  Therefore,  to  third 
order  in  we  have 


This  is  exactly  the  Lioyd-Sholl  result,  and  is  exactly 
the  same  as  the  third  order  energy  based  on  a  linear 
screening  approximation  of  the  self-consistent  pseudo¬ 
potential.  Our  calculation  demonstrates,  however, 
that  it  is  self-consistent  to  third  order  in  the  pseudo¬ 
potential. 
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The  kohn  Vosko  formula  for  electric  field  gradients  associated  with  point  defects  is  modified  to  account  for  the  effect 
of  anisotropic  screenmg  charge  density  on  the  core  enhancement  factor. 


Within  the  past  several  years,  very  powerful  tech¬ 
niques  have  been  developed  that  allow  one  to  deter¬ 
mine.  from  nuclear  magnetic  resonance  experiments 
with  dilute  alloys,  the  magnitude  of  the  electric  Field 
gradient  (EFG)  at  each  of  the  first  few  neighbor  shells 
around  a  solute  ion  [I  5|.  These  EFG  data,  because 
they  relate  directly  to  the  screening  charge  distribu¬ 
tion  induced  by  the  solute  ion.  provide  an  important 
test  of  theories  of  defect  screening  in  metals. 

The  relationship  between  the  screening  charge  den¬ 
sity  and  the  EFG  at  host  ion  sites  was  explored  several 
years  ago  by  Kohn  and  Vosko  [6|  and  Blandin  and 
Friedel  [7],  who  independently  derived  an  approximate 
expression  for  the  EFG.  Their  theory  is  a  first-order 
perturbation  calculation  of  the  screening  charge  den¬ 
sity  and  involves  certain  additional  approximations 
based  on  the  assumption  that  the  impurity -host  ion 
distance  is  large.  Comparison  of  several  calculations, 
based  on  various  versions  of  this  theory,  with  experi¬ 
mental  EFG  data  reveal  certain  consistent  points  of  dis¬ 
agreement  [8-1 1 ) .  The  best  known  of  these  is  that 
calculated  EFG’s  are  almost  always  too  small,  particu¬ 
larly  at  the  nearest-neighbor  shell  where  the  experi¬ 
mental  values  are  usually  larger  by  a  factor  of  two  or 
more.  Tlr  “  .  generally  been  attributed  to  the  ex¬ 
pected  failure  of  the  theory  at  short  distances  [12]. 
Another  possibility  is  that  additional  contributions 
to  the  EFG  from  nonvanishing  net  charge  densities  in 
neighbouring  ion  cells  may  be  appreciable  [13|. 

Recently,  however,  a  new  point  of  disagreement 


has  emerged  that  is  more  difficult  to  explain.  Ac¬ 
cording  to  the  experimental  results  for  aluminum  re¬ 
viewed  by  Minier  and  Ho  Dung  [10],  the  EFG  at  the 
second  neighbor  shell  is  consistently  much  smaller  than 
that  at  the  first,  third,  and  in  some  cases,  even  the 
fourth  and  fifth  shells  for  a  large  number  of  different 
impurity  ions.  As  Holtham  and  Jena  [12]  point  out. 
this  depression  of  the  second-neighbour  EFG  must  be 
a  property  of  the  host  crystal,  and  not  the  solute  ion, 
because  it  happens  in  the  same  way  for  solutes  with 
positive,  negative  and  vanishing  valence  differences. 
However,  their  attempt  to  explain  this  trend  by  refine¬ 
ment  of  the  First -order  theory  was  not  successful. 

This  failure  of  the  first-order  theory  led  Minier  and 
Ho  Dung  [10]  to  suggest  that  the  screening  charge  den¬ 
sity  in  aluminum  is  anisotropic.  If  this  is  the  case,  then 
higher-order  terms  in  the  host  crystal  potential  are 
needed,  because  all  other  terms  in  the  perturbation 
expansion  lead  to  an  isotropic  density.  In  other  words, 
one  should  use  Bloch  functions  as  the  unperturbed 
states,  rather  than  plane  waves,  as  is  done  in  the  first- 
order  theory. 

In  this  communication  we  examine  one  of  the  con¬ 
sequences  of  extending  the  theory  of  the  EFG  beyond 
a  first-order,  plane-wave  treatment.  We  will  show  that 
when  this  is  done,  the  formal  expression  for  the  EFG 
is  similar  to  the  Kohn-Vosko  result,  but  with  a  gener¬ 
alized  definition  of  the  core  enhancement  factor  that 
appears  in  their  formula  [6], 

Our  starting  point  is  the  general  expression  for  the 
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Table  1 

Core  enhancement  correction  factors  for  various  host  ion 
directions. 


Direction  [110|  (100j  {21 1 J  [111] 

1.10  0.80  0.92  0.98 


EFG  at  a  host  ion  site  Rn  in  terms  of  the  screening 
charge  density  5p( x). 


</(*«)=  -2/ 


5p(Rn  +  r')P2(t i  ) 


d-V  , 


(1) 


where  p'  is  the  cosine  of  the  angle  between  and  r , 
and  P2(p  )  is  a  Legendre  polynomial.  This  can  be  written 
in  terms  of  an  approximate  charge  density  8ps(x' )  as 
follows: 


rbp(Rn  +  r)P2(p) 

<l(Rn)= -2a  {R„.kF)f-  ,3  - dV,{2) 


where  a  is  a  ratio  of  integrals  defined  by  eqs.  (1)  and 
(2).  At  this  point,  if  no  approximations  are  introduced, 
the  choice  of  8ps  obviously  has  no  effect  on  the  final 
result,  although  it  does  determine  the  consistent  defi¬ 
nitions  of  the  core  enhancement  factor  a  and  the 
approximate  screening  charge  density.  Now,  however, 
following  Kohn  and  Vosko,  we  will  introduce  an  ap¬ 
proximation  that  eliminates  the  dependence  of  a  on 
host  ion  distance  Rn,  and  another  approximation  in 
the  evaluation  of  the  remaining  integral  in  eq.  (2). 

Because  of  these  approximations  the  final  result  will 
depend  on  our  choice  of  5ps. 

We  choose  to  let  5ps  be  the  smooth  part  of  the 
screening  charge  density  as  defined  in  the  theory  of 
pseudopotentials  [14).  If  we  let  4>®{x)  be  the  pseudo 
wave  function  for  the  perfect  crystal,  0t(x)  the  perturbed 
pseudo  wave  function,  and  i Z'*(x),  0fc(jr)  the  corre¬ 
sponding  true  wave  functions,  then,  to  first  order  in 
the  defect  pseudopotential  W,  the  true  screening  charge 
density  is 


6p(x)  =  2 


k 2  -  k ’2 


*M+c-c-], 


while  the  smooth  part,  6 p.,  is  given  by  a  similar  expres¬ 
sion  with  0®  in  place  of  0*.  Next,  following  Kohn  and 
Vosko,  we  write  the  Bloch  functions  0?  and  0?  in  the 
form 


ip°(x)  =  e‘k'*ut(x) , 


where  uk  is  periodic,  and  evaluate  6 p  and  6 p,  to  lowest 
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6  P(x)' - 

(2n)4|xl- 


order  in  i/|x|.  The  result  is 

nkcQ}  n  n 

— ^  l  k'*  F  F  ^  *  F.-  (•*)  +  c  .c . ) , 

where  f2  is  the  crystal  volume,  ki:x  =  kFx/\x\,  and 

Again,  the  expression  for  6ps  is  similar.  Substitution  in 
the  expression  for  a  leads  to  the  approximate  result 

,  !\^r’)\2P2(p)A\'lr'^ 

a(kr)  =  —  -  --  - .  (3) 

Sl<t>^F(r')]2P2(p  )dir'lr'i 

To  evaluate  the  remaining  integral  in  eq.  (2)  we  ex¬ 
pand  5 ps(x)  in  plane  waves  and  integrate  to  obtain 


/ 5 ps(Rn  +  r  )  P 2(p  )d  V//3 


=  - 1 rrSfi ps(q)P2(q -Rn) e*  ~  -  4  n8ps(Rn ) . 

Q 

where  q  and  Rn  are  unit  vectors.  In  the  last  step  we 
assume  that,  because  Rn  is  large,  the  only  significant 
values  of  q  are  those  parallel  and  antiparallel  to  Rn,  in 
which  case  P2(q  •  Rn)  =  1 . 

The  final  expression  for  the  EFG  is  therefore  similar 
to  the  Kohn-Vosko  formula 

=  I  ™(*f)6ps(*„)  -  (4) 

where  a  is  given  by  eq.  (3)  and  6 ps  is  the  smooth  part 
of  the  screening  charge  density  as  determined  to  first 
order  in  W  using  0®  as  the  unperturbed  wave  functions. 
It  is  easily  seen  that  if  one  uses  the  plane-wave  approxi¬ 
mation  to  0^  and  the  asymptotic  approximation  to  5 os. 
then  eqs.  (3)  and  (4)  reduce  to  the  Kohn-Vosko  for¬ 
mula. 

To  investigate  the  numerical  significance  of  our 
modification  of  the  core  enhancement  factor  we  used 
the  aluminum  form  factor  data  of  Bertoni  et  al.  [15], 
in  a  first-order  perturbation  approximation  to  0*,  to 
calculate  the  ratio  a  /a,  where  a  is  the  enhancement 
factor  in  the  plane-wave  approximation.  The  resulting 
expression  is 

a/a=  [l  +  2  £  — - -J»,(co**£)|_I, 

L  A  #0*2  -|*F-Ar|2 

where  AT  is  a  reciprocal  lattice  vector,  is  the  host 
crystal  form  factor,  and  0^  is  the  angle  between 
and  2AF  -  K.  Numerical  results  are  given  in  table  1 . 
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These  data  suggest  that  the  modification  of  a  re¬ 
quired  by  a  Bloch  wave  treatment  of  defect  scattering 
can  be  significant.  In  the  particular  approximation  used 
here  the  principal  effects  are  a  depression  of  the  en¬ 
hancement  factor  in  the  [100]  direction  and  a  slight 
enhancement  in  the  [1 10]  direction.  These  effects  are, 
however,  much  too  small  to  account  for  all  of  the  ex¬ 
perimentally  observed  depression  of  the  EFG  at  the 
second  neighbor  shell  in  aluminum.  Hopefully,  this 
will  be  explained  by  a  second-order  calculation  of  the 
anisotropic  screening  charge  density,  the  second  factor 
in  eq.  (4),  which  is  now  in  progress. 

This  work  was  supported  by  the  U.S.  Army  Research 
Office. 
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ABSTRACT 

Numerical  computations  of  self-consistent,  second-order  screening  charge 
densities  are  presented  for  a  proton,  an  aluminum  ion,  and  substitutional  zinc 
and  magnesium  in  jellium  with  the  electron  density  of  aluminum.  Comparisons 
with  more  accurate  density  functional  calculations  show  that  a  second-order 
treatment  fails,  as  expected,  for  the  proton,  but  gives  reasonable  results  for 
the  aluminum  ion.  Second-order  corrections  to  the  linear  screening  charge 
densities  are  found  to  be  small  for  the  weaker-scattering  substitutional  de¬ 
fects.  It  is  concluded  that  a  second-order  theory  is  adequate  for  substi¬ 
tutional  defects  with  small  valence  differences,  and  may  find  limited  use  even 
for  strong-scattering  defects,  such  as  the  vacancy,  provided  high  accuracy  is 


not  required 


I .  INTRODUCTION 


The  success  of  the  linear  screening  approximation  in  the  theory  of 
perfect,  non-transition  metals  crystals  has  been  attributed  to  the  small 
values  of  the  pseudopotential  or  model  potential  matrix  elements  that  enter 
such  calculations1.  For  a  perfect  crystal,  matrix  elements  are  evaluated 
only  at  the  rather  large  wavenumbers  corresponding  to  reciprocal  lattice 
vectors,  where  their  values  are  indeed  usually  much  less  than  the  Fermi  energy 
of  the  metal.  This  justifies  the  use  of  first-order  perturbation  theory  in 
the  calculation  of  screening  charge  densities  which,  in  its  self-consistent 
form,  comprises  the  linear  screening  approximation. 

This  argument  cannot  be  applied  to  an  imperfect  crystal,  because  a 
defect  screening  calculation  involves  matrix  elements  evaluated  at  small  wave- 
numnbers,  where  values  can  be  of  the  order  of  the  Fermi  energy.  Hie  linear 
screening  approximation  is  therefore  of  questionable  validity  in  such  prob¬ 
lems,  and  one  must  look  to  a  higher  order  theory.  In  general,  this  presents  a 
formidable  mathematical  difficulty  because  the  problem  of  self-consistency  is 
then  non-linear. 

In  recent  years,  this  nonlinear  screening  problem  for  point  defects  has 
been  the  subject  of  several  numerical  studies.  Most  such  calculations  have 
been  based  on  density  functional  theory^-?,  though  a  few  investigations 
using  the  Korringa-Kohn-Rostoker  ( KKR)  Green* 6-function  method8-10  and  the 
supercell  model1 1-1^  have  also  been  reported.  While  each  of  these  methods 
can  claim  high  accuracy  in  the  treatment  of  nonlinear  screening  effects,  all 
suffer  from  computational  complexities  which  tend  to  restrict  the  class  of 
defect  problems  that  can  be  treated.  In  particular,  it  seems  to  be  imprac¬ 
tical,  at  the  present  time,  to  apply  such  methods  to  anything  more  complex 
than  the  isolated  point  defect,  or,  in  the  case  of  the  supercell  method,  a 
periodic  array  of  point  defects. 
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For  certain  strong-scattering  cases,  such  as  the  screening  of  inter¬ 
stitial  hydrogen,  a  fully  self-consistent  numerical  treatment  is  probably 
essential  for  most  purposes.  On  the  other  hand,  there  are  many  other  defect 
problems  where  the  perturbing  potential  is  not  so  strong  and  where  an  approxi¬ 
mate  treatment  might  suffice,  These  include  substitutional  point  defects, 
stacking  faults  and  similar  planar  defects,  and  perhaps  even  vacancies  and 
interstitials  in  applications  where  high  accuracy  is  not  required. 

We  are  led,  therefore,  to  consider  simpler,  alternate  methods  for  ob¬ 
taining  approximate  solutions  to  the  nonlinear  screening  problem.  In  partic¬ 
ular,  in  this  article  we  examine  the  use  of  second-order  perturbation  theory 
in  the  calculation  of  screening  charge  densities.  Though  the  theory  is,  in 
principle,  applicable  to  an  arbitrary  arrangement  of  scattering  centers,  our 
numerical  results  are  limited  to  point  defects  in  the  jellium  model.  Our  aim 
is  to  provide  some  insight  as  to  the  magnitude  of  second-order  terms  in  both 
strong-  and  weak-scattering  cases,  the  role  of  self-consistency,  and,  in  gen¬ 
eral,  the  validity  of  an  approximate,  second-order  theory  of  defect  screening. 

The  next  section  is  essentially  a  restatement  of  the  second-order 
screening  theory  of  Lloyd  and  Sholl^  and  Brovman,  et  al.,1^  in  terms  of 
defect  and  host  crystal  potentials.  It  serves  mostly  as  background  for  the 
discussion  of  numerical  results,  which  are  presented  in  Sec.  III.  The  calcu¬ 
lations  include  approximate  and  self-consistent  screening  charge  densities  for 
a  proton,  an  isolated  aluminum  ion  (to  simulate  the  vacancy  potential)  and 
substitutional  magnesium  and  zinc.  We  conclude  that  an  approximate  second- 
order  theory,  in  which  the  defect  potentials  are  screened  to  first  order,  is 
adequate  for  weak-scattering,  substitutional  defects,  and  may  prove  useful 
even  for  stronger  scatterers  such  as  vacancies. 
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II.  THEORY 

The  development  that  follows  is  based  on  the  model  potential  theory  of 
screening  by  valence  electrons,  in  which  the  ionic  model  potential  is 
nonlocal1®.  However,  following  Bertoni,  et  al.,1^  we  assume  at  the  outset 
that  the  nonlocal  potential  can  be’  replaced  by  the  local  potential  that  leads 
to  the  same  screening  charge  density  as  a  fully  nonlocal  calculation  in  the 
linear  screening  approximation.  The  Fourier  transform  of  the  smooth  part  of 
the  electron  density,  i.e.,  that  part  derived  from  the  model  wave  function,  is 
then1* 

2  ^  ^ 

Ps<^}  =  8^  1-Gull  W(q)  +  3Ig(pi'P2'q)W(*l>W(*2)  (1) 

Pl*0 

where  p2  =  -(pi  +  q) ,  e  (q)  is  the  dielectric  function,  G(q)  is  the  exchange- 
correlation  correction1,  g(p1,p2*‘l)  is>  the  function  defined  by  Lloyd  and 
sholl1*,  and  w(q)  is  the  Fourier  transform  of  the  self-consistent  potential 
in  Rydbergs. 

The  electron  potential  energy  due  to  valence  electrons  is  related  to  the 
electron  density  through  Poisson's  equation,  modified  to  account  for  exchange 
and  correlation  energies.  Thus  we  have1 

w<q>-wion<<l>  "  p  <1-G(q))(ps(q)  +  pd(q))  (2) 

where  W^on(q)  is  the  ionic  model  potential,  and  p  (q)  is  the  transform  of  the 


depletion  charge  density 
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To  derive  the  second  order  self-consistency  condition  on  W(q)  we 
eliminate  p&(q)  from  Bqs.  (1)  and  (2)  to  obtain 

W(q)  *  WA{q)  +  ^7!  l9<P1*P2»<*)W<P1)W(p2)  (3) 

p^O 

where  W^(q)  is  the  linear  screening  approximation  to  W(q),  i.e., 

W^q)  =*  Wion(q)/e(q)  (4) 

with 

Wion(q)  “  wion<q>  +  ^7  |_1-<3Cq)  ]pd(q)  (5) 

Finally,  from  Eqs.  (1)  and  (3)  we  rewrite  the  density  equation  in  the  more 
convenient  form, 

0,(1)  -  t.  1~c<ql  f  ».(q)  *  W(ql  ~  Vql1  (6) 

8n  1-G(q)  L  1-e(q)  J 

in  which  the  first  term  is  the  linear  screening  approximation  and  the  second 
term  is  the  second-order,  nonlinear  correction. 

Equations  (3)  and  (6)  comprise  a  statement  of  the  general  second-order, 
self-consistent  screening  problem  for  local  potentials.  An  approximate 
solution,  valid  to  third  order  in  Wj^(q},  is  obtained  by  a  single  iteration  of 
Bq.  (3),  i.e.  by  substituting  W^(q)  in  the  sum  on  the  right  side18*  This 
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leads  to  a  relatively  simple,  closed  form  approximation  to  the  screening 
charge  density  which  can  be  applied  to  an  arbitrary  arrangement  of  ions.  The 
accuracy  of  this  approximate  treatment  of  self-consistant  screening  is  one  of 
the  points  we  study  here. 

In  calculations  of  the  defect  screening  charge  density  we  are  concerned 
with  the  change  in  electron  density  caused  by  the  creation  of  a  defect  in  a 
host  crystal.  We  therefore  take 

W(q)  =*  Wh(q)  +  Wd(q)  (7) 

h  d  ■* 

where  W  (q)  is  given  by  Eq.  (3)  for  the  host  crystal,  and  W  (q)  is  the  change 

in  potential  due  to  the  defect.  Similar  expressions  hold  for  the  first  order 

potential  W^(q)  and  for  the  bare  ion  potential. 

Substitution  of  Eq.  (7)  in  Eq.  (6)  gives,  for  the  change  in  ps(q)  due  to 
the  defect 

6p  <£)  =  ^  1-£(q)  f  W  d(|)  +  Wd(q)  -  w/(q)  1 

"  8n  1-G(q)  L  ^  1-E(q)  J  (8) 


(9) 
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h  4  <3  + 

The  last  term  in  Eq.  (9),  which  contains  W  (q)  as  well  as  W  (q),  demonstrates 

the  dependence  of  the  self-consistent  defect  potential,  and  thus  the  screening 
charge  density,  on  the  arrangement  of  neighboring  ions* 

Because  the  host  crystal  potential  is  weak,  the  accuracy  of  this  approx¬ 
imate  theory  is  determined  largely  by  the  relative  magnitudes  of  the  first  and 
second  terms  on  the  right  side  of  Eqs.  (8)  and  (9).  In  general,  one  would 
expect  the  theory  is  useful  only  if  the  second-order  term  in  Eq.  (8)  is  small 
compared  to  the  linear  screening  term,  as  this  may  be  taken  to  imply  that 
third-  and  higher-order  terms  are  smaller  still,  and  can  therefore  be  ne¬ 
glected.  For  the  computations  described  in  the  next  section,  we  therefore 
drop  the  host  crystal  term  in  Eq.  (9)  and  examine  the  changes  in  screening 
charge  density  caused  by  second-order  defect  terms.  This  amounts  to  the 
jellium  approximation,  in  which  case  Eqs.  (8)  and  (9)  reduce  to 


6ps(q)  = 


2 

a _ 

8  tiN 


1-e(q) 

1-G(q) 


<*>A(q)  + 


to(q)  -  (^(q) 
1-  e(q) 


] 


(10) 


and 


arfqj-w^q)  +  ~  Jgfp, /P2/9 )  **>(p1 )  o)(p2>dp1  (11) 

where  w(q)  is  the  defect  form  factor,  Qq  is  the  average  volume  per  ion  and  N 
is  the  number  of  ions  in  the  crystal. 
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III.  NUMERICAL  RESULTS 

Computations  were  performed  for  a  proton,  an  isolated  aluminum  ion  using 
the  Rasolt-Taylor^ ,  Shaw1®  and  Heine-Abarenkov19  model  potentials,  and  for 
substitutional  zinc  and  magnesium  in  an  aluminum  host  using  the  Shaw 
potential.  In  all  cases,  the  averaged  local  potential  defined  by  Bertoni, 
et  al1^,  was  used  in  computing  the  linear  screening  form  factor  wj^(q).  For 
substitutional  zinc  and  magnesium  Wj^(q)  was  taken  as  the  difference  between 
the  impurity  and  host  form  factors,  using  the  impurity  model  potential  param¬ 
eters  and  depletion  charges  calculated  by  Taut  and  Paasch^O  for  an  aluminum 
host  crystal.  Unless  otherwise  noted,  the  exchange  and  correlation  function 
G(q)  was  taken  from  Toigo  and  Woodruff^1. 

The  nonlinear  integral  equation  for  the  self-consistent  form  factor, 

Eq.  (11),  was  solved  by  successive  approximations  using  u)j^(q)  as  the  initial 
approximation.  Computations  were  performed  in  a  71  point  mesh  with  spacing 
Ax  =  0.05  where  x  =  q/2kf,  kj  being  the  Fermi  wavenumber.  Convergence  of  the 
iterative  process  is  illustrated  in  Figure  1  where  we  show  the  linear  screen¬ 
ing  approximation,  which  includes  the  depletion  charge  term  in  Eq.  (5),  and 
approximations  to  w(q)  after  one,  two  and  three  iterations.  For  the  weaker 
difference  potentials  used  for  substitutional  defects  only  one  iteration  was 
required,  while  convergence  of  the  process  for  the  proton  potential  required 
five  iterations. 

In  all  cases,  the  principal  effect  of  nonlinear  screening  is  the  addi¬ 
tion  of  a  repulsive  potential  at  small  values  of  q.  This  causes  a  decrease  in 
the  absolute  magnitude  of  w(q)  for  the  proton  and  aluminum  ion,  as  shown  in 
Figure  1,  and  a  slight  increase  in  the  repulsive  difference  potentials  for 


zinc  and  magnesium.  The  difference  between  the  self-consistent  form  factor 


e 


and  the  form  factor  obtained  after  one  iteration  was  found  to  be  small  for 
aluminum,  somewhat  greater  for  the  proton,  and  negligible  for  substitutional 
zinc  and  magnesium. 

The  fact  that  the  form  factor  is  altered  only  at  smaller  values  of  q 
deserves  comment.  From  Eq.  (10)  we  see  that  the  second-order  correction  term 
is  divided  by  1-e(q)  which  diverges  in  the  long  wavelength  limit.  The  corre¬ 
sponding  limit  of  6ps(q),  which  is  proportional  to  the  integrated  screening 
charge,  is  therefore  unaffected  by  the  addition  of  the  second  order  term,  as 
is  required  by  charge  conservation.  Also,  because  the  difference  u>(q)  -a)j^(q) 
is  nonvanishing  only  at  small  q,  where  e(q)  is  large,  Eq.  (10)  tells  us  that 
the  Fourier  transform  of  the  screening  charge  density  is  not  as  sensitive  to 
nonlinear  effects  as  is  the  form  factor  itself.  lhis  may  explain  the  obser¬ 
vation  that  nonlinear  corrections  to  the  screening  charge  density  are  less 
pronounced  than  corrections  to  other  functions  derived  from  the  form  factor, 
such  as  the  interionic  potential  and  phonon  spectra^.  It  should  also  serve 
as  a  warning  that  the  conclusions  drawn  from  comparisons  of  charge  densities 
may  not  apply  to  the  calculation  of  other  properties. 

Calculated  screening  charge  densities  are  shown  in  Figures  2  through  6. 
For  the  proton  and  aluminum  ion  we  show  the  linear  screening  result,  the 
approximate  second  order  density  and  the  self-consistent  second  order  density. 
For  substitutional  zinc  and  magnesium  the  differences  between  the  approximate 
and  self-consistent  distributions  are  negligible. 

The  results  for  the  proton,  which  are  presented  in  Figure  2,  show  a 
strong  nonlinear  effect  and  a  significant  difference  between  the  approximate 
and  self-consistent  solutions.  At  large  distances  the  principal  effects  are  a 
large  shift  in  phase  and  an  increase  in  the  amplitude  of  the  Friedel  oscilla¬ 
tions,  while  at  small  distances  there  is  a  pronounced  increase  in  electron 
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y.  Although  these  trends  are  in  qualitative  agreement  with  the  results 
e  accurate  calculations,  quantitative  agreement  is  not  obtained, 
ations  based  on  the  density  functional  approach^  show  a  somewhat 
r  shift  in  phase  and  a  smaller  amplitude  in  the  long  range  oscillations, 
greater  increase  in  density  near  the  proton.  At  the  proton  site,  for 
e,  the  value  shown  in  Figrure  2  is  about  one-third  the  correct  value. 
Results  obtained  with  the  Shaw  potential  for  an  aluminum  ion  are  shown 
[ure  3.  Here  again  we  find  that  second  order  effects  are  significant, 
lat  changes  in  the  screening  charge  density  are  not  as  dramatic  as  was 
ise  for  a  proton.  The  principal  effects  are  a  slight  shift  in  the  phase 

1  increase  in  amplitude  of  the  Friedel  oscillations  with  some  distortion 

2  density  inside  the  core.  In  this  case,  the  difference  between  the 
consistent  and  approximate  treatments  of  second  order  screening  is 
3ible  at  most  distances. 

We  also  performed  self-consistent  calculations  using  the 
t-Taylor^  and  Heine-Abarenkov^  potentials  for  the  aluminum  ion.  The 
ts  are  presented  in  Figure  4,  which  also  includes  the  self-consistent 
lata.  Another  calculation  based  on  the  Rasolt-Taylor  potential  with  the 
ctric  function  used  by  Dagens,  et  al*  is  not  shown  here  because  the 
t  does  not  differ  significantly  from  the  Rasolt-Taylor  curve  in  Figure  4. 
nclude,  from  a  comparison  of  Figures  3  and  4,  that,  except  at  short  dis- 
s,  the  changes  produced  by  second-order  screening  are  of  roughly  the  same 
tude  as  those  associated  with  the  use  of  different  model  potentials. 

It  is  perhaps  more  useful  to  compare  our  second-order  approximation  for 
Luminum  ion  with  the  more  accurate  density  functional  calculations  ra¬ 
il  by  Manninen,  et  al^.  From  their  Figure  1  and  our  Figures  3  and  4 
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we  see  that/  at  distances  greater  than  2.0  a.u.,  our  results  lie  roughly  in 
between  their  curves  for  the  jellium  model  and  the  jelliura  vacancy  model 
(jellium  with  the  positive  background  removed  from  the  ion  cell)/  though  the 
amplitude  of  our  Friedel  oscillations  is  slightly  greater.  While  this  does 
not  prove  that  a  second-order  treatment  is  adequate,  it  does  at  least  show 
that  screening  charge  densities  calculated  in  this  approximation  are  reason¬ 
able,  even  for  strong  scattering  centers  like  isolated  ions  and  vacancies. 

Finally,  in  Figures  5  and  6  we  show  the  first-  and  second-order  charge 
densities  for  two  weak-scattering  centers,  substitutional  zinc  and  magnesium. 
Here  again  the  principal  nonlinear  effects  are  a  shift  in  phase  and  an  in¬ 
crease  in  amplitude  of  the  Friedel  oscillations,  and  a  rearrangement  of  charge 
in  the  core.  These  effects  are,  however,  quite  small  compared  to  the  non¬ 
linear  effects  obtained  for  the  proton  and  aluminum  ion,  presumably  because 
the  potentials  used  here  were  derived  from  differences  in  ion  potentials  and 
are  therefore  weaker  than  potentials  for  isolated  ions.  If  this  is  so,  then 
calculations  for  substitutional  impurities  with  greater  valence  differences 
should  show  stronger  nonlinear  effects. 

IV.  CONCLUSION 

Self-consistent,  second-order  screening  charge  densities  have  been  cal¬ 
culated  in  the  jellium  approximation  for  a  proton,  an  isolated  aluminum  ion 
and  substitutional  magnesium  and  zinc  in  aluminum.  As  one  would  expect,  the 
results  do  not  agree  with  more  accurate  density  functional  calculations  for 
the  strongest  scattering  center,  the  proton,  although  the  changes  produced  by 


the  addition  of  second-order  terms  are  qualitatively  correct.  For  the  substi¬ 
tutional  defects,  numerical  results  indicate  that  a  simpler  theory,  in  which 
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form  factors  are  calculated  in  the  linear  screening  approximation,  is  suitable 
provided  the  valence  difference  is  small.  In  the  case  of  somewhat  stronger 
perturbations,  such  as  the  isolated  ion  or  vacancy,  the  accuracy  of  a  second- 
order  treatment  is  subject  to  question.  Still,  the  changes  in  the  screening 
charge  distribution  produced  by  the  self-consistent  addition  of  second-order 
terms  are  reasonable,  and  the  theory  might  therefore  find  limited  application 
even  in  such  strong  scattering  situations. 
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FIGURE  CAPTIONS 


Figure  1 


Figure  2 


Figure  3 


Figure  4 


Figure  5 


Figure  6 


.  Self-consistent  form  factor  for  aluminum.  'Ate  lower  solid  curve  is 
the  linear  screening  approximation  including  depletion  charge 
effects,  while  the  upper  solid  curve  is  the  second-order  self- 
consistent  form  factor.  The  dotted  and  dashed  curves  were  obtained 
after  one  and  two  iterations,  respectively,  of  Eq.  11. 


Second-order  screening  charge  density  for  the  proton.  The  dotted 
curve  is  the  linear  screening  approximation,  the  dashed  curve  cor¬ 
responds  to  the  approximate  second-order  form  factor  obtained  after 
one  iteration  and  the  solid  curve  is  the  self-consistent  result. 


Second-order  screening  charge  density  for  the  aluminum  ion.  The 
legend  is  the  same  as  in  Figure  2. 


Self-consistent,  second-order  screening  charge  densities  for  the 
aluminum  ion.  The  dotted  curve  corresponds  to  a  damped 
Heine-Abarenkov  potential1®,  the  dashed  curve  to  the  Rasolt-Taylor 
potential3  and  the  solid  curve  to  the  optimized  model  potential16. 


Screening  charge  density  for  substitutional  magnesium  in  an 
aluminum  host.  The  dotted  curve  is  the  linear  screening  result 
while  the  solid  curve  is  the  self-consistent,  second-order  density. 


Screening  charge  density  for  substitutional  zinc  in  an  aluminum 
host.  The  legend  is  the  same  as  in  Figure  5. 
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FIGURE  4 
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ABSTRACT 

Electric  field  gradient  (EFG)  calculations  for  metal  alloys  are  usually 
based  on  an  asymptotic  theory  in  which  the  valence  effect  EFG  is  proportional 
to  the  screening  charge  density.  We  show,  by  means  of  second  order  computa¬ 
tions  that  avoid  the  usual  asymptotic  approximation,  that  this  proportionality 
does  not  hold  for  hetrovalent  alloys  of  aluminum.  In  such  cases,  the  charge 
density  is  oscillatory  and  is  significantly  affected  by  the  addition  of  aniso¬ 
tropic  second-order  terms,  while  the  EFG  is  dominated  by  the  isotropic  first 
order  term  and  is  non-oscillatory .  Hiis  failure  of  the  usual  asymptotic  EFG 
theory  is  attributed  to  the  inadequacy  of  the  asymptotic  approximation  in  the 
long  wavelength  limit.  An  approximate  correction  term  is  derived  and  found  to 
be  proportional  to  the  difference  in  valence  between  impurity  and  host  ions. 
Additional  first-order  computations  show  that  the  corrected  EFG  formula  gives 
good  agreement  with  results  obtained  by  direct  numerical  evaluation  of  the 
exact  expression  for  the  EFG. 


I. 


INTRODUCTION 


When  a  point  defect  is  created  at  a  host  ion  site  in  a  cubic  metal,  the 
resulting  disturbance  produces  electric  field  gradients  ( EFG ) ,  which  are 
measurable  by  nuclear  magnetic  resonance  methods, at  neighboring  host  ion 
sites.  As  is  well  known,  the  EFG  has  two  components,  the  so-called 
valence-effect  EFG,  caused  by  the  redistribution  of  conduction  electrons,  and 
the  size-effect  EFG,  caused  by  lattice  distortion. 3'4  The  valence-ef feet  EFG 
is  generally  assumed  to  be  given  by  the  Kohn-Vosko  formula^  (also  derived 
independently  by  Blandin  and  Friedel®)  which  states  that  the  EFG  at  a  given 
host  ion  site  is  directly  proportional  to  the  screening  charge  density  at  that 
site.  Calculations  of  the  size-effect  component  are  based  on  the  theory  of 
Sagalyn,  et  al,3  which  assumes  that  the  EFG  tensor  is  proportional  to  the 
local  strain  tensor. 

As  was  demonstrated  by  Sagalyn  and  Alexander4  and  others,7-^  both  the 
valence-  and  size-effect  contributions  must  be  included  to  obtain  satisfactory 
agreement  with  experimental  EFG  data.  Realizing  this,  we  nevertheless  chose 
to  concentrate  on  the  valence-effect  EFG  in  the  work  reported  here.  The 
reason  we  made  this  choice  is  that  our  original  intent  was  to  investigate 
anisotropic  screening  effects  in  the  hope  of  explaining  the  anomalously  small 
values  of  the  EFG  observed  at  the  second  neighbor  site  in  most  aluminum 
alloys. ^0  This  we  sought  to  do  through  a  nonlinear,  second-order  calculation 
of  the  screening  charge  density,  as  explained  in  Section  II. 

The  result  of  our  calculation,  which  is  also  presented  in  Section  II, 
shows  a  pronounced  anisotropy  in  the  screening  charge  distribution.  In 
Section  III,  however,  we  show  that  if  we  avoid  an  asymptotic  approximation  in 
the  evaluation  of  the  valence-effect  EFG  integral,  then  second-order  screening 
terms,  which  include  the  anisotropic  part,  are  insignificant  compared  to  the 


rder,  linear  screening  term  in  cases  where  the  impurity  and  host 
s  differ.  As  a  corollary  to  this  result  we  find,  contrary  to  the 
■sko  formula,  that  the  EFG  is  not  proportional  to  the  screening  charge 
’  in  such  hetrovalent  alloys. 

lis  leads  to  the  main  point  of  the  article,  which  we  present  in  Section 
;re  we  reexamine  the  approximations  involved  in  the  Kohn-Vosko  treatment 
rive  a  correction  term,  proportional  to  the  valence  difference,  to  their 
Stic  formula.  Through  numerical  computations,  we  then  show  that  the 
ted  asymptotic  formula  gives  good  agreement  with  results  obtained  by 
numerical  integration.  We  conclude,  therefore,  that  when  host  and 
ty  valences  differ,  the  correction  term  derived  here  is  essential  to  the 
te  prediction  of  the  valence-effect  EFG  in  simple  metal  alloys. 

II.  ANISOTROPIC  SCREENING  CALCULATIONS 
n  the  preceding  article11  we  reported  numerical  studies  of  the 
onsistent,  second-order  screening  charge  density  for  various  point 
s  in  aluminum.  The  calculations,  which  were  performed  in  the  jellium 
imation,  indicate  that  a  second-order  treatment  should  suffice  for 
cattering  substitutional  defects  for  which  the  difference  between  the 
nd  solute  valences  is  not  too  large.  In  this  section  we  extend  the 
ations  to  include  second-order  scattering  events  involving  the  host 
1. 

s  before,  we  define  the  defect  potential  as  the  difference  in  model 
ials  between  the  impurity  and  host  ions.  When  the  model  potential  for 
st  crystal  is  added  to  the  jellium  model,  the  screening  charge  density 
s  anisotropic  and  the  self-consistent  defect  potential  is  therefore  also 
ropic.  However,  based  on  recent  self-consistent  supercell  calculations 
vacancy  in  aluminum, ^  the  anisotropy  in  the  defect  potential  is  expected 
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to  be  weak.  We  therefore  average  the  potential  over  direction  and  thus  obtain 
the  following  integral  equation  for  the  self-consistent  defect  form  factor 
c*Xq ) : 

Uq)  =  ub(q)  +  “£  -^/g(q,K,p)w(p)dQp  (1) 

ic*o 


with 

uo(q)  =  u)ji(q)  +  f^L1  ]^Z9(q'P'r)  w<P)a^r>  (2) 

q  q  + 

P 

In  these  expressions  c(q)  is  the  Hartree  dielectric  function,  G(q)  is  the 
exchange-correlation  correction,  is  the  host  crystal  form  factor  at 

reciprocal  lattice  vector  K,  g(q,K,p)  is  the  function  defined  by  Lloyd  and 
Sholl,13  and  wj^(q)  is  the  defect  form  factor  in  the  linear  screening 
approximation.  The  integral  in  Eq.  (1)  is  over  the  direction  of  the  vector 
p  =  -(q+k),  and  r  =  |q+p|  in  Eq.  (2). 

The  Fourier  transform  of  the  smooth  part  of  the  screening  charge  density 
is 

&Ps(q)  =  6p°(q)  +  g(q,K,p)u£w(p)  (3) 

K+0 


where 


6p®(q) 


oMqJ-w  (q) 
1-e(q) 
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(4) 


The  isotropic  part  of  the  screening  charge  density,  6ps(q),  has  the  same  form 
as  in  the  jellium  model.  However,  it  does  not  equal  the  density  obtained  in 
the  jellium  approximation  because  the  potential  u^(q)  that  appears  in  Eq.  (4) 
is  coupled  to  the  spherical  average  host  crystal  potential  through  Eqs.  (1) 
and  (2). 

For  numerical  computation  of  the  screening  charge  density,  we  derived 
local  form  factors  from  the  exponentially-damped,  nonlocal  potentials  of  Heine 
and  AbarenV.ov^ ^  by  the  averaging  method  of  Bertoni,  et  al,^5  and  used  data 
from  Toigo  and  Woodruff *6  for  the  exchange-correlation  function  G(q).  The 
self-consistency  condition,  Eq.  (1),  was  solved  by  an  iterative  method 
described  elsewhere,^  and  the  resulting  self-consistent  defect  potential  was 
then  used  to  compute  the  transform  of  the  screening  charge  density  according 
to  Eqs.  (3)  and  (4).  The  inverse  transform  is  shown  in  Figure  1  for  two 
crystallographic  directions,  along  with  the  isotropic  term  6p°. 

These  results  indicate  that  the  addition  of  the  host  crystal  potential  to 
the  jellium  model  •  coduces  a  significant,  anisotropic  distortion  of  the 
screening  charge  at  large  distances  from  the  defect,  and,  as  shown  in  the 
inset  in  Figure  1 ,  an  isotropic  shift  in  the  smooth  part  of  the  charge  density 
at  the  impurity  site.  This  rather  large  change  at  the  impurity  site  is 
consistent  with  the  calculations  of  Manninen  and  Nierainen,17  who  found  that 
the  addition  of  a  spherical  average  crystal  potential  to  their  density 
functional  model  of  a  vacancy  also  produced  a  large  shift  in  density  at  the 
vacancy  site.  Our  principal  interest  at  present,  however,  is  in  the 
anisotropic  screening  charge  density  at  neighboring  ion  sites.  If  indeed  the 
EFG  at  these  sites  is  proportional  to  the  smooth  part  of  the  screening  charge 
distribution,  then  the  results  displayed  in  Figure  1  say  that  the  addition  of 
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the  host  crystal  potential  should  produce  anisotropic  shifts  in  the  EFG 
distribution  over  neighboring  ion  shells.  Ibis  is  the  point  we  address  in  the 
next  section. 

III.  SECOND-ORDER  CONTRIBUTIONS  TO  THE  EFG 

If  we  let  6p(r)  be  the  total  change  in  valence  electron  density,  which 
includes  changes  in  the  oscillatory  part  inside  host  ion  cores,  then  the 
valence-effect  EFG  at  host  ion  position  R  is  given  by4»5,9 

qv(  R)  =  -  +  ~ne6p( R)  +  2e/P2(r*R)  6p(R+r)^j^- 


=  “  +  |-ne6p(R)  -  |jp£6p(q)P2(q*R)e 

-*• 

q 


(5) 


where  r,  R  and  q  are  unit  vectors,  P2  (cos©)  is  a  Legendre  polynomial,  and  AZ 
is  the  impurity  valence  minus  the  host  valence.  This  is  the  principal 

component  of  the  traceless  valence  effect  EFG  tensor  as  defined  by  Sagalyn  and 

.  V 

Alexander;4  in  their  notation  the  quantity  we  calculate  is  — eq | |  .  Our 

analysis  differs  from  that  of  Ponnambalam  and  Jena^  because  our  charge  density 

is  not  spherically  symmetric  with  respect  to  the  impurity  ion,  which  requires 

that  we  use  the  general  form  given  by  Eq.  (5),  rather  than  the  spherically 

symmetric  form  derived  by  Ponnambalam  and  Jena. 

Next,  following  Kohn  and  Vosko,^  we  assume  that  6p(q)  can  be  replaced 
by  a6pg(q)  where  a  is  the  core  enhancement  factor  that  provides  an  approxi¬ 
mate  correction  for  the  effects  of  charge  density  oscillations  in  the  core  and 
6pg  is  the  smooth  part  of  the  valence  electron  density.  From  the  original 


6 


I 

development,^  It  can  be  seen  that  the  core  enhancement  approximation  rests  on 

|  the  assumption  that  R  is  large,  i.e.  it  is  an  asymptotic  approximation.  We 

■* 

I  might,  therefore,  also  assume  that  the  only  important  q  directions  in  the 

integral  in  Eq.  (5)  are  those  that  are  nearly  parallel  to  R,  in  which  case 
!  P2(q.R)  ~  1.  In  this  case  the  integral  becomes  proportional  to  6ps(R) 

and  we  recover  the  Kohn-Vosko  formula.  (This  will  be  demonstrated  in  a 
slightly  different  way  in  the  next  section.)  However,  in  the  calculations 
reported  below  we  chose  to  avoid  this  particular  approximation  because  it  was 

A  A 

a  rather  simple  matter  to  insert  the  factor  P2(q*R)  in  the  numerical 

integration  routine  already  developed  for  our  charge  density  calculations, 

and  thus  provide  a  test  of  at  least  one  aspect  of  the  large  R  approximation. 

Typical  results  are  shown  in  Table  I.  The  entries  are  EFG  contributions 

in  units  of  10 13  cgs-esu  for  a  core  enhancement  factor  a  -  10.0.  The  major 

column  headings.  Term  1,  etc.,  correspond  to  the  first,  second  and  third  terms 

in  Eq.  (5).  The  subheadings  1,  2  and  3  correspond,  respectively,  to  the 

linear  screening,  second-order  isotropic,  and  second-order  anisotropic 

■* 

components  of  the  charge  density  6ps(q). 

The  most  important  point  to  be  noted  here  is  that  the  integral  term 
(Term  3)  is  not  proportional  to  6ps(R).  In  fact,  the  linear  screening 
contribution  to  this  term  t Column  1),  which  dominates  at  all  R,  is 
non-oscillatory,  while  6ps(R)  itself  exhibits  the  familiar  Friedel 
oscillations  at  large  R.  The  immediate  conclusions  are  that  (1)  the 

A  A 

approximation  p2(q.R)~1  leads  to  serious  error,  and  (2)  nonlinear  screening 
to  the  EFG  are  not  as  important  as  one  might  expect  from  calculations  of  the 


charge  density  alone 


We  are  led,  therefore,  to  a  reexamination  of  the  asymptotic  approximation 


as  it  affects  the  integral  term  in  Eg.  (5),  and  also  as  it  affects  the 
core  enhancement  approximation,  which  is  also  based  on  the  assumption  that  R 
is  large.  In  the  development  that  follows  we  will  drop  nonlinear  terms 
because,  as  can  be  seen  from  Table  I,  the  linear  screening  part  of  the 
integral  term  is  the  dominant  component  of  the  EFG  at  all  R. 

IV.  A  NEW  ASYMPTOTIC  APPROXIMATION 


We  return  to  Eq.  (5)  and  focus,  for  the  moment,  on  the  total  change  in 
charge  density,  6p(r).  To  avoid  the  core  enhancement  approximation  we  must 
determine  the  impurity-induced  change  in  the  total  valence  electron  wave 
function,  including  the  oscillatory  part  in  the  cores  of  host  ions.  To  do 
this  we  use  linear  screening  theory  and  the  core  state  projection  operator18 
in  a  straightforward,  first-order  perturbation  calculation.  The  resulting 
expression  for  the  charge  density  is 

+  -*• 


**•*  -5^7 


-G(q) 


(6) 


q 


where  f(q,r),  which  determines  the  shift  in  orthogonalization  charge  density, 
is  given  by 


q2i1-e<q) 

8nL1-G(q) 


]f(q,r)  =  l 


f  (k) 
o 


t  *2-lit*;l2 


[<k|  1-p|rXr|  1-p|k+q> 


-  <k|rXr | k+q>  J 


(7) 
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where  fQ(k)  is  the  electron  occupation  number  for  electron  state  |)c>,  and  p 
is  the  core  state  projection  operator  for  a  host  ion  at  the  origin.  In  the 
Appendix  we  show  that  substitution  of  Eq.  (6)  in  the  integral  term  in  Eq.  (5) 
leads  to 

4eQ  00  2 

qint(R)  =  ~Yr!  Iril-I-G(q)  ^(q^1+F(q)  ^j2(qR)q2<3q  (8) 

o  H 

where  C^j  is  the  volume  per  ion,  j2<qR)  is  a  spherical  Bessel  function,  and 
F(q)  is  the  function  defined  in  the  Appendix  and  plotted  Figure  2. 

a  a 

The  approximation  P2(q.R)~1  in  Eq.  (5)  is  equivalent  to  using  the 
asymptotic  approximation 

j2(qR) - j0<qR>  (9) 

in  Eq.  (8).  If  we  do  this,  and,  in  addition,  replace  F(q)  by  its  value  at  the 
singular  point  q  =  2kf,  then  Eq.  (8)  reduces  to  the  Kohn-Vosko  form 

8  TUB 

<Iint(R) - — oc6ps(R)  (10) 

with  a  =  1+F(2kj).  As  we  show  in  the  Appendix,  the  second  term  in  Eq.  (5)  can 
also  be  written  in  this  form.  Thus,  by  combining  the  second  and  third  terms 
in  Eq.  (5),  we  find  that  the  electronic  part  of  the  EFG  is  given  by  Eq.  (10) 
with  a  =  1+F(2kf)-p,  where  J3  is  the  quantity  defined  in  the  Appendix.  The 
numerical  value  of  a  based  on  Herman-Ski liman  1®  core  functions  for  aluminum 
and  the  on-the-Fermi-sphere  (OFS)  approximation^  to  nonlocal  terms  is  6.8, 
which  compares  favorably  with  other  calculations  of  the  core  enhancement 
factor  for  aluminum. 21 


r 

i 

i 

The  difficulty  with  this  development  is  that  Eq.  (9)  becomes  a  poor 
'  approximation  as  q  >0/  and  leads  to  a  serious  inaccuracy  in  the  estimation  of 

qfntOO*  Ibis  is  illustrated  in  Figure  3,  which  is  a  plot  of  the  function 

*(q)  =  j2(qR) 

=  q26ps(q)  j2<qR)  (ID 

and  the  approximation  to  <Mq)  obtained  from  Eq.  (9).  In  this  illustration,  R 
is  the  nearest  neighbor  distance  in  aluminum,  and  u(q)  is  the  average  defect 
form  factor  for  a  magnesium  impurity  in  aluminum.  We  note  that,  because  6ps(q) 
is  the  screening  charge  density  for  a  hetrovalent  impurity,  charge  conservation 
requires  that  6ps(q)  be  nonvanishing  as  q  -*■  0.2°  For  a  homovalent  impurity, 
however,  6ps(q)  ■»  0  as  q  +0,  and  the  error  introduced  by  the  asymptotic 
approximation  is  not  as  large  as  that  shown  in  Figure  3.  Ibis,  we  shall  see, 
is  reflected  in  the  correction  term  derived  below. 

To  obtain  a  better  approximation  to  qj.nt(R)  we  add  the  term 
4eQ  “  2  i_  /  l 

£q(R)  =  “2^  8^U^G(q)  ]w(q)[1+F(q)]ji(qR)qdq  (12) 

where  we  have  used  the  identity  jQ(Z)  +  j2<Z)  =  3ji(Z)/Z.  A  very  simple  form 
results  if  we  assume  that  the  integrand  in  this  expression  is  negligible  for 
large  q,  take  the  long  wavelength  limits  of  q2e(q)  and  F(q),  and  cut  the 
integration  off  at  qg,  the  first  zero  of  u>(q).  With  M(q)  defined  by 

u><q)  =  -  TT^l  M(q)  03) 

wrfl 

we  obtain 

£q(R)  ~  4eAZ^'t'F(?)J  J^qM(q) j ^ (qR)dq  (14) 

o 


I 
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We  now  assume  that  £q(r)  is  insensitive  to  the  form  of  M(q)  as  long  as 
M(0)  *  1  and  M(qg)  *  0.  To  obtain  an  integrable  form  we  therefore  take 

M(q)  ~j0(7q/q0>  (15) 

which  gives,  to  lowest  order  in  1/R, 

Aj(R)  ~  2eAZL1+F(°)  J/r3  (16) 

A  new  asymptotic  formula  for  the  EFG  is  obtained  by  combining  &i(R)  with  our 
estimates  of  the  other  terms  in  Eq.  (5).  The  result  is 

q(R)  ~  2e^f  -  (°-  -  “[l+F(2kF)-p]6ps(R)  (17) 

Figure  4  shows  the  two  components  of  q  as  well  as  the  total  EFG  according 
to  Eq.  (17),  along  with  data  obtained  by  direct  numerical  integration  of 
Eq.  (5),  again  for  magnesium  in  aluminum.  Similar  results  were  obtained  for 
substitutional  zinc,  silicon  and  germanium  in  aluminum.  For  gallium,  with 
AZ=0,  the  familiar  oscillatory  form  was  obtained  in  both  calculations  with 
agreement  between  the  two  being  somewhat  poorer  than  that  shown  here.  Thus, 
it  appears  that  the  asymptotic  approximation  leading  to  the  first  term  in 
Eq.  (17)  is  adequate,  and  that  the  result  is  reasonably  accurate  when  AZtO  and 
the  first  term  dominates. 


11 


V.  CONCLUSION 

Equation  (17)  is  the  principal  result  of  the  present  work.  It  is  a 
corrected  asymptotic  formula  for  the  valence  effect  EFG  caused  by  an  impurity 
ion  with  valence  difference  AZ.  The  correction  to  the  Kohn-Vosko  formula,  the 
first  term  on  the  right  side  of  Eq.  (17),  is  the  sum  of  the  EFG  due  to  the 
bare  impurity  ion  and  the  EFG  associated  with  the  long  wavelength  behavior  of 
the  screening  charge  density.  For  large  impurity-host  ion  distance,  which  is 
assumed  here,  the  ionic  term  cancels  exactly  with  the  plane  wave  part  of  the 
screening  term.  This  leaves,  as  the  only  contribution  from  the  long 
wavelength  term,  the  EFG  caused  by  the  shift  in  crthogonalization  charge 
associated  with  a  completely  screened  impurity  ion.  The  result  is 
proportional  to  AZ  and  is  the  dominant  contribution  to  the  EFG  unless  the 
valence  difference  vanishes. 

The  development  leading  to  Eq.  (17)  shows  why  nonlinear  screening  does 
not  contribute  significantly  to  the  valence  effect  EFG  even  though  nonlinear 
effects  are  important  in  calculations  of  the  charge  density  itself.  This  is 
because  when  Az+O  the  dominant  contribution  to  the  EFG  comes  from  the  long 
wavelength  limit  of  the  charge  density  or,  equivalently,  from  an  integral  of 
the  charge  density  over  all  space,  and,  as  is  well  known,  this  limit  is  given 
exactly  by  linear  screening  theory. 20  m  other  words,  nonlinear  corrections 
vanish  as  q  +  0  and  therefore  do  not  affect  the  first  term  in  Eq.  (17). 

Although  our  numerical  results  pertain  only  to  alloys  of  aluminum,  we 
believe  the  principal  conclusion  to  be  valid  for  all  simple  metal  alloys.  The 
heart  of  the  matter-  is  the  failure  of  the  usual  asymptotic  approximation  in 
the  long  wavelength  limit,  as  illustrated  in  Figure  3.  Regardless  of  how  one 
corrects  for  this  deficiency,  the  result  must  reflect  the  behavior  of  the 
charge  density  at  small  values  of  the  wavenumber  and  must,  therefore,  be  at 
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APPENDIX 

To  calculate  f(q,r)  we  expand  the  plane  wave  and  core  functions  in 

Eq.  (7)  in  spherical  harmonics  in  a  coordinate  system  with  the  polar  axis  in 

the  q  direction.  Then,  with  the  help  of  the  addition  theorem  for  spherical 

harmonics  and  the  recursion  relation  for  Legendre  polynomials  to  simplify 

terms  involving  2p  core  states,  we  find  that  f(q,r)  is  a  function  of  the 

magnitudes  of  q  and  r  and  \l,  the  cosine  of  the  angle  between  q  and  r*  This 

-*  •* 

allows  us  to  write  f(q,r)  as  a  sum  of  Legendre  polynomials  in  |i.  Then,  after 
substitution  in  Eqs.  (6)  and  (5)  and  integration  over  r  we  obtain  the 
following  result  for  the  integral  term  in  Eq.  (5): 

-►  ■* 

8ne  (  q± 

qmt(R)  -  -  ^rJ  ^ 

q 

a  a 

where  q  and  R  are  unit  vectors  and  F(q)  is  given  by 

F(q)  =  -  f/°°  7  J1  f(q#r,n)P  ((i)d^idr 
o  1  -1  * 

If  we  use  the  on-the-Fermi-spheri  (OFS)  approximation^  to  the  matrix 

■* 

elements  in  Eq.  (7)  then  the  k  integral  is  easily  evaluated  and  the  resulting 
expression  for  F(q)  is 

F(q)  ■  )^I(2L+1 )  [fuL<kF,q)  +  ffli/hq/q)  ]  +  7 
NL 


L7if{f}]w<q)Vq‘R,e  l1+F(q)J 
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where 


*«.<*'*>  -  St-n^NLOOn  ir"*<2*HM2r+1)(*  *  r)2P«(5q) 

u' 


rV<r\..,  .  aL  . 

*L  —p—U  <qr)d(icr7L:,-ltkrldr 


2.,-  P2.<r’ 


For  q  <  2kf,  kq  =  |k+q|  =  kf  and  £q  =  q/2kf,  while  for  q  >  2kf,  kq  =  q  -  kf 
and  =  1  (OFS  approximation).  The  numbers  N  and  L  are  core  state  quantum 

numbers  and  PflL(r)  is  the  corresponding  radial  wave  function.  To  calculate 
the  function  shown  in  Figure  2,  we  used  the  core  functions  tabulated  by  Herman 

and  Skillman.^9 

The  quantity  p  that  appears  in  Eq.  (17)  is  determined  by  the  function 
f(q,r)  evaluated  r  =  0.  From  Bq.  (7)  it  is  easily  shown  that,  in  the  OFS 
approximation 

f<q,o)  .  -  (*F>  *  ('„„<*,>  -  r„«vj  Vo'WoKo 


where 


R„  =  iim  P  (r)/r 
NL  r^o  NL 

If  we  approximate  f(q,0)  by  its  value  at  q  =  2kf  then  we  find,  from  Eq.  (6), 


6p(R)  ~  2 p6ps(R) 


where 


P  -  +  f (2k  ,0)  ] 
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Terns  1 ,  2  and  3  refer  to  the  first,  second  and  third  terms  in  Eq.  (5),  while  subheadings  1,  2,  and  3  correspond 
respectively,  to  linear  screening,  second-order  isotropic  and  second-order  anisotropic  components  of  the  charge 
density*  The  letter  a  denotes  the  lattice  constant. 
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FIGURE  CAPTIONS 


Figure  1.  Anisotropic  screening  charge  density  in  the  110  and  112  directions 

for  a  magnesium  ion  in  aluminum.  Hie  isotropic  term  6ps°  is  defined 
by  Eq.  (4). 

Figure  2.  Hie  function  F(q)  that  governs  the  orthogonalization  charge 
contribution  to  the  electric  field  gradient  as  in  Eq.  (8). 

Figure  3.  The  solid  curve  is  the  function  <|>(q)  defined  byEq.  (11);  the  dashed 
curve  is  the  approximation  obtained  by  using  Eq.  (9).  Hie  extra 
node  in  the  approximate  curve  leads  to  a  large  error  in  the  electric 
field  gradient  integral  unless  solute  and  solvent  valences  are 
equal. 

Figure  4.  Contributions  to  the  electric  field  gradient  for  a  magnesium 

impurity  in  aluminum.  Hie  lower  dotted  curve  is  the  first  term  in 
Eq.  (17),  the  upper  dotted  curve  is  the  second  term  and  the  dashed 
curve  is  their  sum.  Hie  solid  curve  is  the  result  obtained  by 
numerical  integration  of  Eq.  (8). 
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APPENDIX  H 


Third-Order  Calculation  of  Lattice 
Distortion  for  a  Vacancy  in  Aluminum 


(Draft  of  a  paper  to  be  submitted  to  Physical  Review  B) 
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In  this  report  we  examine  the  influence  of  third-order  forces  on  the 
lattice  distortion  field  for  a  vacancy  in  aluminum.  He  chose  to  do  this 
particular  calculation  because  our  earlier  study  of  second-order  corrections 
to  the  screening  charge  density,  which  give  rise  to  third  order  forces, 
produced  results  in  reasonable  agreement  with  more  accurate  density  functional 
calculations.  There  was,  therefore,  some  reason  to  expect  that  a  third-order 
force  approximation  would  suffice,  at  least  as  an  indicator  of  the  relative 
magnitudes  of  second-  and  third-order  terms. 

For  our  calculation  we  used  the  method  of  lattice  statics  with  the 
second-  and  third-order  force  expression  given  by  Solt  and  Werner  for  a 
hetrovalent  alloy  with  potential  difference  Av.  For  the  vacancy  we  took  Av 
to  be  the  negative  of  the  aluminum  ion  potential  in  the  linear  screening 
approximation.  The  potential  used  for  this  purpose,  and  for  calculating  the 
dynamical  matrix,  was  that  developed  by  Sen,  et  al.,  with  screening  given  by 
Taylor's  dielectric  function.  This  potential  is  a  two-parameter  model 
potential  of  the  Heine-Abarenkov  form,  with  parameters  adjusted  so  as  to 
obtain  the  best  overall  agreement  with  a  number  of  experimentally  determined 
properties  of  aluminum,  including  the  phonon  spectrum.  Because  the  potential 
of  Sen,  et  al,  gives  good  agreement  with  the  phonon  spectrum  in  a  second-order 
calculation,  we  did  not  include  third-order  terms  in  our  calculation  of  the 
dynamical  matrix  and  its  inverse,  the  lattice  Green's  function. 

With  the  lattice  statics  method,  the  calculation  of  lattice  displacement 
vectors  involves  an  integration  over  the  first  Brillouin  sons.  For  this 
purpose  we  used  the  special-direction  method  with  a  64  point  gaussian 
integration  along  each  of  21  directions  in  1/48th  of  the  sone.  In 
second-order,  the  displacements  obtained  from  this  integration  differed  by 
less  than  2%  from  those  obtained  with  32  point  integrals  in  10  directions. 
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Second-  and  third-order  displacements,  in  units  of  one-half  the  lattice 
constant,  are  given  in  Table  I.  Our  second-order  displacements  are  much 
smaller  than  those  reported  by  Singhal,  particularly  in  the  first  four 
neighbor  shells.  This  we  attributed  to  the  use  of  different  pseudopotentials. 
The  addition  of  third-order  forces  makes  most  of  the  displacements  smaller 
still,  and  even  produces  sign  changes  for  the  second  and  fifth  neighbors. 

This  indicates  that  third-order  ion-vacancy  forces  are  generally  repulsive, 
tending  to  decrease  the  inward  relaxation  of  ions  around  the  vacancy,  in 
agreement  with  Hasegawa's  finding  that  third-order  interionic  forces  are 
generally  attractive. 

The  disturbing  point  illustrated  by  these  data  is  the  magnitude  of 
third-order  effects,  particularly  at  the  first  and  second  neighbor  sites.  If 
we  are  to  believe  that  third-order  terms  are  actually  as  large  as  indicated 
here,  then  we  must  be  concerned  about  the  magnitudes  of  fourth-  and  higher 
order  terms  that  were  not  included  in  our  calculation.  In  other  words,  we 
must  then  question  the  usefulness  of  a  third-order  theory  of  interionic  forces 
involving  defects.  On  the  other  hand,  our  calculation  is  a  local 
approximation,  and  this  could  have  led  to  significant  error  in  the  calculation 
of  forces.  We  tend  to  believe  that  the  second-order  part  of  the  calculation, 
including  the  determination  of  the  dynamical  matrix,  is  reasonably  accurate 
because  it  is  based  on  an  empirical  potential  fitted  to  experimental  data 
through  second-order  theory.  However,  the  success  of  this  local  potential  in 
second-order  calculations  does  not  justify  its  use  in  third-order  terms  where 
nonlocal  effects  may  be  quite  different. 

Had  third-order  effects  been  smaller,  we  might  be  justified  in  ignoring 


higher-order  terms  and  the  uncertainties  introduced  through  a  local 
approximation  for  the  third-order  term.  It  is  clear,  however,  that  this 


cannot  be  done.  We  oust  conclude,  therefore,  that  the  local,  third-order 
theory  of  interionic  forces  involving  defects  is  subject  to  serious  question. 
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